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The  subject  of  the  present  effort  is  numerical  developments  in  support  of  a  BiGlobal  linear  insta¬ 
bility  analysis  [26]  of  the  compressible  three-dimensional  steady  laminar  basic  flows  obtained  within 
the  framework  of  F61775-99-WE049.  The  first  objective  becomes  extraction  from  the  aforementioned 
results  of  the  appropriate  two-dimensional  basic  flow  data.  The  assumption  is  made,  and  is  verified 
in  the  numerical  results  to  be  increasingly  applicable  as  the  flow  Mach  number  increases,  that  the 
dependence  of  the  basic  flow  on  an  appropriately  defined  downstream  spatial  direction  is  much  weaker 
than  that  in  the  wall-normal  and  the  circumferential  directions  of  the  elliptic  cone.  In  terms  of  the 
basic  flow,  derivatives  in  the  downstream  direction  are  thus  neglected;  in  the  disturbance  equations  an 
eigenmode  Ansatz  is  introduced  in  this  direction.  A  discussion  is  presented  regarding  the  choice  of  coor¬ 
dinate  system  for  the  analysis.  The  viscous  disturbance  equations  of  the  two-dimensional  compressible 
partial-derivative  eigenvalue  problem  are  presented  in  cartesian  coordinates.  It  is  subsequently  pro¬ 
posed  that  on  grounds  of  numerical  efficiency  the  analysis  be  performed  in  an  inviscid  framework  on 
the  elliptic  confocal  coordinate  system.  The  governing  equation  is  derived  in  the  chosen  coordinate 
system  and  numerical  means  for  its  discretisation  are  discussed. 

In  an  inviscid  global  linear  analysis  the  issue  of  avoidance  of  singularities  of  the  governing  equations 
through  indentation  of  the  path  of  integration  must  be  addressed.  However,  the  concept  of  direction 
is  ill-defined  in  the  framework  of  an  inviscid  global  analysis.  One  alternative  to  shooting  in  both  the 
local  and  the  global  theory  is  to  pose  the  disturbance  equations  as  a  matrix  eigenvalue  problem  on 
complex  calculation  grids.  A  new  spectral  collocation  method  has  been  devised,  which  employs  this 
technique  and  is  validated  by  solving  several  one-dimensional  linear  eigenvalue  problems  on  different 
families  of  collocation  grids.  Tensor  products  of  the  optimal  grids  identified  through  the  present  effort 
can  be  employed  for  the  solution  of  the  two-dimensional  eigenvalue  problem  at  hand,  if  marginally 
amplified/damped  or  neutral  modes  are  sought.  The  suite  of  numerical  tools  proposed  for  the  BiGlobal 
instability  analysis  has  been  validated.  First,  model  Poisson  problems  have  been  solved  on  the  elliptic 
confocal  coordinate  system  and  the  expected  exponential  convergence  property  has  been  demonstrated. 
Subsequently,  three  flow  configurations  in  the  subsonic,  supersonic  and  hypersonic  regimes  have  been 
solved.  Most  interestingly,  large-scale  non-azimuthally-periodic  disturbances,  inaccessible  to  classic 
linear  theory,  have  been  identified  in  all  flow  regimes. 
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1.  Introduction 

The  present  effort  concentrates  on  global  linear  instability  analyses  of  the  steady  laminar  flows 
on  an  elliptic  cone  presented  in  [25].  These  flows  are  intrinsically  three-dimensional;  as  a  matter  of 
fact  the  three-dimensionality  of  the  geometry  of  the  elliptic  cone  accounts  for  the  departure  of  the 
laminar-turbulent  transition  process  from  the  fairly  well- understood  scenaria  on  circular  cones  [21] 
and  other  axisymmetric  bodies  of  revolution.  Shedding  light  on  the  essentially  three-dimensional  tran¬ 
sition  process  on  the  elliptic  cone  from  the  perspective  of  global  linear  analysis  provides  motivation 
for  several  current  experimental  [15,  16,  12,  6]  and  numerical  [29,  30]  efforts,  as  well  as  for  the  theo¬ 
retical/numerical  methodology  described  herein. 

The  most  general  framework  in  which  a  global  linear  instability  analysis  can  be  performed  is  that 
of  three-dimensional  linear  theory  in  which  all  three  spatial  directions  are  resolved.  This  is  consistent 
with  the  separability  in  the  governing  equations  of  time  on  one  hand  and  the  three  spatial  directions 
on  the  other.  Formulation  of  the  three-dimensional  global  linear  eigenvalue  problem  is  straightforward; 
however,  its  numerical  solution  is  not  feasible  with  present-generation  computer  architectures.  Indeed, 
coupled  resolution  of  dim  spatial  directions  requires  storage  of  two  arrays  for  the  real  and  imaginary 
parts  of  the  eigenfunctions,  each  of  size 

5  x  26  x  N2xdim  x  10“9  Gbytes 

of  core  memory  in  primitive  variables  formulation  and  64-bit  arithmetic  if  N  points  resolve  each 
spatial  direction.  The  size  of  each  array  is  doubled  if  128-bit  arithmetic  is  deemed  to  be  necessary 
[7].  If  a  numerical  method  of  optimal  resolution  power  for  a  given  number  of  discretisation  points  is 
utilised,  such  as  spectral  collocation,  experience  with  the  one-dimensional  eigenvalue  problem  suggests 
that  in  excess  of  N  —  64  must  be  used  for  adequate  resolution  of  eigenfunction  features  at  critical 
Reynolds  numbers  typical  of  boundary  layer  flows.  The  resulting  sizes  of  the  respective  matrices  are 

~  22  Tbytes  ~  5.4  Gbytes  and  ~  1.3  Mbytes, 

if  dim  —  3,  2  and  1,  respectively.  It  becomes  clear  that  while  the  classic  linear  local  analysis  ( dim  — 
1)  requires  very  modest  computing  effort  and  is  indeed  part  of  industrial  prediction  toolkits,  the 
main  memory  required  for  a  three-dimensional  linear  instability  analysis  is  well  beyond  any  currently 
available  or  forecast  computing  technology.  Between  the  two  limits,  a  global  linear  analysis  based 
on  solution  of  the  two-dimensional  linear  eigenvalue  problem  is  well  feasible  and  has  indeed  been 
performed  using  a  variety  of  numerical  methods  and  platforms,  delivering  results  inaccessible  to  local 
theory  (e.g.  [22,  24])  or  new  insights  into  old  instability  problems  (e.g.  [28]);  a  review  is  presented  in 
[27]. 

It  thus  becomes  clear  that  in  order  to  proceed  with  a  global  analysis  of  the  elliptic  cone  basic  flows 
the  disturbance  equations  must  be  posed  in  the  framework  of  solution  of  a  two-dimensional  eigenvalue 
problem.  The  underlying  assumption  is  that  the  basic  flow  is  independent  of  one  spatial  direction;  in 
the  problem  at  hand  an  appropriate  direction  must  be  identified  such  that  the  dependence  of  the  basic 
flow  on  it  is  much  weaker  than  that  on  the  other  two  spatial  directions.  While  the  arc-length  direction 
is  an  obvious  candidate  to  be  considered  as  homogeneous  in  the  global  instability  problem,  the  choice 
of  the  other  two  spatial  directions  is  not  unique.  Indeed,  three  alternative  directions  complementing 
the  arc-length  and  wall-normal  coordinates  are  identified  in  what  follows.  From  the  point  of  view  of 
numerical  feasibility,  the  elliptic  confocal  coordinate  system  is  one  natural  basis  on  which  to  express 
and  collocate  the  various  forms  of  the  two-dimensional  compressible  eigenvalue  problems. 

The  viscous  two-dimensional  eigenvalue  problem  is  derived  herein  in  cartesian  coordinates.  From  it 
the  generalised  Rayleigh  equation  presented  in  [25]  may  be  derived;  the  latter  is  expressed  here  in  the 
elliptic  confocal  coordinate  system.  If  no  loss  of  physical  information  is  to  be  expected  it  is  attractive 
to  consider  working  with  the  matrices  discretising  the  inviscid  global  eigenvalue  problem,  whose  size 
is  (3/5)2  that  of  their  viscous  counterparts.  However,  a  new  issue  is  raised  when  working  in  an  invis¬ 
cid  framework,  namely  the  appearance  of  critical  layers  and  potential  singularities  of  the  governing 
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equations.  This  problem  is  well-known  in  local  linear  theory,  where  its  solution  through  indentation 
of  the  path  of  integration  is  straightforward.  This  is  so  because  there  is  no  ambiguity  in  defining  the 
path  of  integration;  in  other  words,  the  line  on  which  basic  flow  data  is  defined  in  a  boundary-layer 
context  is  the  wall-normal  spatial  direction.  The  necessary  basic  flow  data  are  obtained  on  the  in¬ 
dented  path  by  Taylor  expansion  around  the  critical  layer.  In  such  a  framework,  efficient  numerical 
approaches  exist,  based  on  finite-differences  and  shooting  [8],  to  obtain  the  leading  eigenvalue.  In  a 
global  analysis  on  the  elliptic  cone,  on  the  other  hand,  the  basic  flow  is  defined  on  a  plane  and  both 
concepts  of  integration-path  and  shooting  are  inapplicable.  This  led  us  to  propose  a  new  algorithm 
for  the  integration  of  the  inviscid  disturbance  equations,  based  on  complex  spectral  collocation  grids. 
Technical  details  are  presented  by  Theofilis  et  al.  [13].  Tensor  products  of  the  optimal  grids  identified 
through  the  present  effort  will  be  employed  for  the  solution  of  the  two-dimensional  eigenvalue  problem 
at  hand. 

Progress  has  been  made  on  several  fronts.  Firstly,  in  §  2  we  present  for  the  first  time  the  viscous 
compressible  partial-derivative  eigenvalue  problem  in  cartesian  coordinates.  The  inviscid  limit  of  the 
eigenvalue  problem  is  then  recalled  [25]  as  being  the  most  efficient  method  of  analysis  from  a  numerical 
point  of  view.  The  elliptic  confocal  coordinate  system  is  thus  introduced  in  §  3,  which  on  the  one  hand 
preserves  the  numerical  advantage  of  efficiency  and  on  the  other  offers  natural  means  of  resolving 
boundary-layer  flow  features  by  an  analytic  exponentially  stretched  body-fitted  coordinate  system. 
Secondly,  we  also  discuss  in  this  section  the  continuing  quest  to  express  the  equations  of  motion 
in  alternative  coordinate  systems.  An  approximation  of  the  type  made  to  derive  the  two-dimensional 
eigenvalue  problem  in  cartesian  coordinates  must  be  made  subsequently  and  it  is  interesting  to  compare 
the  respective  systems  with  respect  to  optimal  representation  of  the  underlying  physical  instability 
mechanisms.  Ultimately,  experience  from  experiment  could  decide  the  optimal  choice  of  coordinate 
system. 

The  main  thrust  of  the  current  effort  has  been  along  the  first  path  and  in  §  4  we  present  the  numerical 
means  by  which  the  generalised  Rayleigh  equation  may  be  collocated  on  the  elliptic  confocal  coordinate 
system  and  solved  as  either  a  temporal  or  a  spatial  two-dimensional  eigenvalue  problem  at  the  same 
level  of  numerical  effort.  However,  in  view  of  the  issue  of  two-dimensional  critical  layers,  particular 
to  the  inviscid  two-dimensional  eigenvalue  problem,  we  proceed  in  §  5  to  elaborate  on  the  concept  of 
a  complex  calculation  grid.  The  latter  permits  performing  inviscid  instability  analyses  by  solving  the 
eigenvalue  problem  in  a  direct  manner.  Several  validation  cases  in  the  limit  of  one-dimensional  linear 
inviscid  equations  are  presented.  Attention  is  then  focussed  in  §  6  on  the  BiGlobal  eigenvalue  problem. 
First,  the  proposed  spectral  collocation  numerical  discretisation  scheme  is  validated  by  solution  of 
Poisson  problems  on  the  elliptic  confocal  coordinate  system.  Subsequently,  two-dimensional  planes  of 
steady  basic  flow  data  on  different  aspect  ratio  elliptic  cones  are  extracted  from  three-dimensional 
solutions  at  M  —  0.5,4  and  8.  No  analysis  from  a  physical  point  of  view  as  such  has  been  performed 
here;  rather,  attention  has  been  focussed  on  numerical  aspects,  such  as  feasibility  of  the  methodology, 
resolution  and  its  influence  on  different  classes  of  eigendisturbances  and  extent  of  the  integration 
domain  and  its  relation  to  the  boundary  conditions  imposed.  It  is  demonstrated  that  application  of 
BiGlobal  analysis  is  relatively  straightforward  in  terms  of  the  resolution  necessary  for  subsonic  flow, 
although  the  assumptions  underlying  the  derivation  of  the  compressible  generalised  Rayleigh  equation 
are  increasingly  applicable  to  supersonic/hypersonic  flow.  Classes  of  BiGlobal  eigendisturbances  that 
can  be  recovered  accurately  in  high-speed  flow  are  identified.  Interestingly,  results  obtained  provide 
indications  that  BiGlobal  instability  theory  can  be  utilised  as  an  alternative  to  classic  linear  theory 
in  investigations  of  vortex  breakdown  over  delta  wing  configurations  in  high-angle  of  attack  and  as  a 
tool  alternative  to  DNS  in  computational  aeroacoustics  studies.  A  short  discussion  in  §  7  summarises 
the  present  effort,  while  we  use  the  Appendix  to  present  technical  details. 
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2.  The  viscous  compressible  two-dimensional  eigenvalue  problem 

In  this  section  we  start  from  first  principles  to  derive  the  equations  describing  the  compressible  two- 
dimensional  eigenvalue  problem.  This  governs  small-amplitude  perturbations  developing  upon  a  steady 
laminar  compressible  basic  state  which  is  inhomogeneous  in  two  spatial  directions.  The  amplitude 
functions  themselves  are  two-dimensional  functions  of  the  resolved  spatial  directions.  We  choose  to 
focus  the  discussion  on  cartesian  coordinates  in  view  of  their  wider  utility  as  compared  with  any  of 
the  body-fitted  elliptic  cone  coordinate  systems  discussed  in  the  previous  section.  In  this  framework 
any  flow  quantity  may  be  decomposed  as 

q  =  q(®,  y)  +£  q(®,  y)  E,  (2.1) 

where 


E  —  exp  \i{fiz  —  tot) 

and  q  =  {().  u,  v,  w.p)1  is  the  vector  of  the  three-dimensional  nonconservative  variables  embed¬ 
ded  in  cartesian  space  Oxyz7  q  =  (p(x,y),u(x,y),v(x,y),w(x,y),p(x,y))T  is  the  vector  of  the  two- 
dimensional  steady  basic  flow  quantities  and  q  =  (p(x,y),u(x,y),v(x,y),w(x,y),p(x,y))T  is  that  of 
the  amplitude  functions  of  the  three-dimensional  disturbances,  which  are  periodic  in  one  spatial  di¬ 
rection  alone,  z.  As  usual,  an  equation  of  state  for  ideal  gases  relates  pressure  p  and  density  p  through 
temperature  T, 


7M2p  =  pT,  (2.2) 

and  a  general  law  relating  viscosity  p  and  temperature  is  adopted.  Upon  decomposition  of  temper¬ 
ature  and  viscosity  into  their  steady  two-dimensional  and  unsteady  disturbance  components  one  also 
has 


T  =  T(x,y)  +eT{x,y )  E, 

p  =  p(x,y)  +  ep{x,y)E.  (2.3) 


Further  quantities  introduced  are  the  viscous  stresses 


and  the  viscous  heat  fluxes 


T** 

=  (v/Re) 

[4/3 ux  -  2/3{vy  +  wz) 

(2.4) 

Tyl' 

=  (v/Re) 

[4/3 vy  -  2/3 (ux  +  wz) 

(2.5) 

II 

05 

[4/3 wz  -  2/3 (ux  +  vy) 

(2.6) 

txy 

=  tYX  = 

{p/Re){uy  +vx) 

(2.7) 

TXZ 

=  TZX  = 

{p/Re){uz+wx) 

(2.8) 

rYZ 

=  tzy  =  (p/Re){vz  +  wy), 

(2.9) 

1 

II 

* 

O* 

^Tx 

(7  -  1  )M2PrRe 

(2.10) 
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_ pTy _ 

(7  —  1  )M'2PrRe 


pTz 

(7  —  1  )M2PrRe 


In  this  context,  the  equations  of  motion  can  be  written  in  conservative  form, 


(2.11) 

(2.12) 


where 


Pt  +  (pu)x  +  (pv)y  +  ( Hz  =  0 
[pu)t  +  (p  +  pu2)x  +  ( puv)y  +  ( puw)z  -  txx  -  txy  -  txz  =  0 
(pv)t  +  ( puv)x  +  {p  +  pv2)y  +  ( pvw)z  -  TXy  -Tyy  -  TYZ  =  0 
(pw)t  +  ( puw)x  +  ( pvw)y  +  (p  +  pw2)z  -TXZ  -TyZ  -  TZZ  =  0 
et  +  [u(e  +  p)\x  +  [v(e  +  p)\y  +  [w(e  +  p)\z 

+Qx  +  Qy  +  Qz 


—  \UTXX  +  VTXY  +  WTXZ 


—  \uTXY  +  VTY  Y  +  WT 1  Z 


iz  ,  rz  ,  zz 

—  UT  +  VT  +  WT 


=  0 


(2.13) 

(2.14) 

(2.15) 

(2.16) 


(2.17) 


e  = 


7 


— -p  +  ~p(r/2  +  v2  +  w2) 


(2.18) 


and  subscripts  denote  partial  differentiation  with  respect  to  an  independent  variable.  The  global 
viscous  linear  analysis  proceeds  by  substitution  of  (2.1)  into  (2.13-2.17)  upon  which  systems  of  equa¬ 
tions  at  0(1),  0(e)  and  0(e2)  are  obtained  in  a  manner  conceptually  identical  with  that  of  classic 
linear  theory  [17].  The  terms  at  0(1)  represent  the  steady  basic  state,  those  at  0(e2)  are  neglected  in 
a  linear  framework,  while  the  system  of  equations  to  be  solved  at  0(e)  reads 
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Disturbance  Continuity 

vpy  +  vpy  +  upx  +  upx  +  p  (— iu ;  +  i  (3  w  +  vy  +  ux)  +  p  (i  (3  w  +  vy  +  ux)  —  0,  (2-19) 

Disturbance  X-Momentum 


u  v  py  +  u  v  py  +  u  v  py  +  px  +  2  u  u  px  +  u2  px 


+p{vuy  +  u  (—i  oo  +  i  (3  w  +  vy  +  2  ux )  j 
+p  (—i  w  u  +  i  (3  uw  +  i  (3  uw  +  v  uy  +  v  uy  +  uvy  +  uvy  +  2uux  +  2u  ux) 

^  I  ^  3  Py  Uy  3  p  nyy  4~  2  i  f3  w  px  2  Vy  px  ~\~  2  vy  px  4  px  ux 

4  Px  ^x  3  Py  {'U'y  4”  ^ x )  3  Py  ^ x  3  i  (3  p  Wx  ft  Vxy  4  p  llxx 

P  ^  3  (3  u  -\-  3  Uyy  i  [3  wx  +  vxy  4 uxx^j  ^  —  0, 


+  - 


(2.20) 


Disturbance  Y- Momentum 


py  +  2  v  v  py  +  v2  py  +  u  v  px  +  u  v  px  +  u  v  px 
+p  ju  (— i  uj  +  i  f3  w  +  2vy  +  ux )  +««XJ 

+p  (—iuiv  +  if3vw  +  if3vw  +  2vvy  +  2vvy  +  vux  +  vux  +  uvx  +  uvx) 

^  ^  4“  2  i  (3  w  Py  4  py  Vy  4  Py  vy  3  i  (3  p  wy  4  p  vyy  3  uy  px 

3  Uy  fix  “l-  2  fly  Ux  4"  2  fiy  ux  3  fix  Ux  3  fix  Vx  f  Uxy  3  fi  Vxx 

-  p  (—3  (32  V  +  i  [3  wy  +  4  vyy  +  uxy  +  3  flxx)  }  =  0,  (2.21) 
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Disturbance  Z- Momentum 

i  (3p  —  iujpw  +  2i(3pww  +  vwpy  +  vwpy  +  vwpy  +  pwvy+pwvy 
+p  vwy  +  pvwy  +  uwpx  +  uwpx  +  uwpx  +  pwux  +  pwux  +  puwx 
+p  (y—ioow  +  i/3w2  +  wvy  +  vwy  +  wux  +  u  w x^j  +puwx 

-5-  \  —  3i  PvJIy  +  2i  P  p,Vy  —  3  fly  Wy  —  3  JlyWy  —  3  jl  Wyy 
O  il6  ^ 

—3i/3u]Xx  +  2i/3fiux  —  3fixwx  —  3jXxwx  —  3fi  wxx 
—  JJ  ^—4  /32  w  +  i  (3  vy  +  3  Wyy  +  i  P  ux  +  3  w  xx)  } 


+ 


Disturbance  Energy 
i  co 


2(7-l) 


p  (u2  +  v2  +  w2J  +2  p  ( 


|2p  +  (7  —  1)  p(u2  +  v2  +  w2'' 


UU  +  V  V  +  w  w) 


)]} 


+ 


—  I  (7  -  1)  u3  px  +  2'ypux  -  pv2ux  +  7 pv2ux  -  2pvvux 


2  (7-I) 


_ A  _  A  _ O _  A  _ O  _  _  _  A  _ 

+27 pvvux  —  pw  Ux  +  'ypw  ux  —  2pwwux 


+27 pwwux  +  2r)pux  —  pv2  ux  +  7 pv2  ux 


-pw2ux  +  7 pw2ux  +  3  (7  -  1)  u2  ( pux  +  pux ) 

u  [2  7  Pa;  +  (7  -  1) 

^3  u2  px  +  v2  px  +  w2  px  +  6  p  u  ux  +  2  p  v  vx  +  2  ~p  w  wx^j  j  j 
+  2  ^_1j{  +  2'!px-2vvpx  +  2'yvvpx-2wwpx 
+2  7  w  w  px  -  v2  px  +  7  v2  px  -  w2  px  +  7  w2  px 


—2  pvvx  +  27  pvvx  —  2  pvvx  +  27  pvvx 
—2pvvx  +  2rypvvx  —  2pwwx  +  2ripwwx 
—  2~pwwx  +  2'y~pwwx  —  2pwwx  +  2'ypwwx'\ 
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- - —  j  (7  —  1)  V3  f)y  +27  pVy  ~  PI?  Vy  +7/3  Vy  ~  2  P  U  U  Vy 


2  (7-I) 


_ _  _  ^  _2 _  ^  _ o _  _  _  ^  ^  _ 

+  27  pUUVy—pW  Vy+'ypw  Vy  —  2pWWVy+2r)pWWVy 


+  2  7  PVy  —  pu 2  Vy  +  7  pu2  Vy  ~  /9  10  2  Vy  +  7  /J®2  Vy 


+  3  (-1  +  7)  V2  ( PVy  +  PVy) 
V  [2  7  py  +  (7  -  1) 

( ?2py  +  3  v2py  +  w2py  +  2  puUy  +  6  pvijy  +  2pwWyj  j  | 
+  ^ — jy  {2^py  —  2ii,upy  +  2'yuupy  —  2wwpy  +  2'ywwpy 

—  U 2  f)y  +  ”/  U2  f)y  —  W2  Py  +  ^  W2  {)y  —  2  pUUy  +  2  f)UUy 


—  2  p-UUy  +27  pUUy  —2  pUUy  +  2  7  /9  U  Uy  ~  2  jjWWy 
+  27  pw  Wy  —2  pWWy  +27  pWWy  —  2pw  Uly  +  2'ypw  Wy^ 

+  2^  ^  ^  |2 7pro  —  2fluuw  +  2^/~puuw  —  2pvvw 
+2”/pvvw  +  (— 1  +7)  pw  (u2  +  v2  +  w2^ 


+27  pw  —  pu2  w  +  7  pu2  w  —  pv2  w 


_ o  A  _  _ o  A  _  _ o  A  'j 

+  7 pv  w  —  3pw  w  +  3~/pw  w> 

”(7-  1  )M2  Pr-Rei^ Tx  +  /ix^'x  +  ^Txx  +  11  ^xx 

Py  ^'.7  3'  Py  Ty  -\-  jlT yy  -\-  P  Tyy  ^ 

(7-l)M2  Pr  Re 

+  — —  l  +  3v  uypx  +  3vuy~px  +  3  vuy  px  —  2i  /3~pw  ux  —  2  pvy  ux 

3  Pe  l 
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—2jAvyux  +  Afiux  +  3if3jIwux  —  2jIVyUx  +  8jAuxux 
+3  (iuyvx  +  3jIUyVx  +  3vJIxvx  +  3v  Jaxvx  +  3  jlvx 
+3jAuyvx  +  3vJIxvx  +  §jAvxvx  +  3wJIxwx  +  3wfAxwx 
+3  Ja  wx  +  3  w  Jlx  wx  +  6  Jawx  wx  +  3  (a  v  uxy  +  3  JIv  uxy 
-\~3  fi  v  uXy  3  Ja  v  vxx  3  [A  v  vxx  3  fA  v  vxx  3  {a  w  wxx 
+3  Jaw  wxx  +  3  Jaw  wxx  +  3  u  JJy  uy  +  3  u  fAy  uy  +  3  (AUy 
+3  uJIyUy  +  ftJJUyUy  —  2i  (3jJw  vy  +  A  fatty  +  3i  [3Jaw  Vy 

+  8jAVyVy+3wJIyWy  +  3wfiyWy+3fiWy+3wjAyWy 
+  6jAWyWy  +  3  flu  Uyy  +  3  JI  U  Uyy  +  3  JI  U  Uyy  +  3  [A  W  Wyy 

+3  Jaw  wyy  +  3  JIw  wyy  —  2  JAvyux  —  2jAvyux  —  2jZvyux 


+3u  jLyvx  +  3u p,y  vx  +  3  jiUy  vx  +  3  jiuy  vx  +  3  u  jiy  vx 
A~  3  fA  Uy  Vx  ~\~  3  jA  U  Uxy  A~  3  fA  U  Uxy  A~  3  fA  U  Uxy  ^ 

+  3^ei  3^^^  ~  2vyjAx  +  AjJxux  +  3i  j3jJwx  -  2jJvxy  +  AjJuxx} 
A~  |  2  i  J3  w  jix  2  Vy  fix  2vy  (ax  +  A  Jax  ux  +  4  fix  ux 

2  i  f3  /a  wx  2  fi  vXy  2  /a  vXy  ~\~  4  jl  uxx  A~  4  fA  uxx  ^ 


-^3i  /3  w  iAy  +  A/AyVy  +  3i  J3  jAWy  A-  An  vyy  —  2  /Ay  ux  —  2  /a  uxy  j 

v  r  _  _  _  _ 

+  3^Re  >■  “  2  iP  ™  Vy  +  4  %  VV  +  4  Uy  Vy  -  2  iJ3  n  Wy  +  4  Ja  Uyy 

A~  A  fA  Vyy  2  jly  UX  2  fA y  UX  2  fA  Uxy  2  fA  Uxy  ^ 
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—ijjfl 
3  Re 


(Vy  +  Ux) 


( VWy  +UWX)  | 


t*P 
3  Re 


fivv  +  Afiww  +  2iwvy  +  2iwvy  —  3  iv 


wv 


—  3iv  wy  +  2iw  ux  +  2iw  ux  —  3iuwx  +  3u  (/3  u  —  i  wx)  }  =  0 
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These  equations  may  be  recast  as  follows.  First,  in  line  with  analogous  practice  of  classic  linear 
theory,  viscosity  is  taken  to  be  a  function  of  temperature  only.  A  law  fj,  —  fj,(T)  is  considered  to  relate 
the  two  quantities,  taking  viscosity  in  general  to  be  a  nonlinear  function  of  temperature.  A  Taylor 
expansion  about  the  basic  state 


p  —  p(T  +  eT) 


H(f)  +  s^f  +  0(s2) 


(2.24) 


may  then  be  used  to  identify  the  basic  state  viscosity  and  relate  that  of  the  disturbance  field  to  the 
disturbance  temperature  T  through 


(2.25) 

Equation  (2.25)  may  be  used  to  eliminate  ft  from  the  system  (2.19-2.23)  in  favour  of  the  disturbance 
temperature.  The  latter  is  related  with  the  respective  disturbance  pressure  and  density  of  the  medium 
through 


7M2p  =  pT  +  Tp ,  (2.26) 

which  may  be  used  to  eliminate  one  of  p.  p.  T  and  cast  the  final  system  to  be  solved  in  terms  of  five 
two-dimensional  disturbance  amplitude  functions,  the  three  velocity  components  and  two  thermody¬ 
namic  variables.  Once  a  solution  of  the  eigenvalue  problem  is  known  (2.26)  may  be  used  to  calculate 
the  eliminated  third  thermodynamic  disturbance  quantity. 

(a)  Challenges  and  potential  simplifications  of  the  compressible  two-dimensional  eigenvalue  problem 

The  challenges  associated  with  a  numerical  solution  of  the  system  (2.19-2.23)  can  be  classified  into 
two  categories.  One  is  related  with  issues  analogous  with  those  of  the  incompressible  two-dimensional 
eigenvalue  problem  [26]  and  stems  from  the  simultaneous  resolution  of  two  spatial  directions,  in  this 
case  the  problem  being  aggravated  by  the  existence  of  an  additional  constitutive  equation  to  be  solved. 
In  addition,  the  system  to  be  solved  becomes  increasingly  stiff  with  increasing  Reynolds  number; 
note  that  the  diagonal  dominance  of  the  matrix  decreases  linearly  with  Re.  This  leads  to  the  need 
for  simplifications,  which  will  be  discussed  in  what  follows.  The  second  class  of  challenges  is  raised 
by  compressibility.  From  a  numerical  point  of  view  mixed  derivatives  of  both  the  basic  flow  and  the 
disturbance  quantities  appear,  which  couple  the  system  stronger  but  are  prone  to  propagating  errors  in 
the  numerical  discretisation  in  one  spatial  direction  into  the  entire  resolved  plane.  Cubic  nonlinearities 
give  rise  to  the  need  for  finer  discretisation  of  the  compressible  eigenvalue  problem  compared  with  its 
incompressible  counterpart.  An  additional  source  of  stiffness  in  the  eigenvalue  problem  stems  from 
the  quadratically  decreasing  diagonal  dominance  of  the  system  of  equations  with  increasingly  high 
supersonic  and  hypersonic  Mach  number.  On  the  positive  side,  the  issue  of  boundary  conditions  for 
the  (elliptic)  incompressible  partial-derivative  eigenvalue  problem  is  absent  in  compressible  flow,  where 
pressure  becomes  a  variable  to  be  solved  for  as  opposed  to  a  constrain  to  be  imposed. 

In  dealing  with  the  issue  of  the  size  of  the  eigenvalue  problem  one  notes  that,  unlike  the  incom¬ 
pressible  eigenvalue  problem  where  considering  a  wavenumber  vector  perpendicular  to  the  plane  on 
which  the  basic  flow  develops,  i.e.  setting  w  =  0,  leads  to  the  ability  to  recast  the  eigenvalue  problem 
as  one  with  real  coefficients  requiring  approximately  half  the  computing  resources  for  its  numerical 
solution  compared  with  the  original  complex  system,  in  compressible  flow  this  simplification  is  not 
possible.  This  is  so  on  account  of  density  and  viscosity  variations  which  are  absent  in  incompressible 
flow.  Focussing  on  the  application  at  hand,  even  if  such  a  simplification  were  possible  it  would  be 
physically  meaningless,  since  the  predominant  flow  direction  on  the  elliptic  cone  coincides  with  that  of 
the  wavenumber  vector.  Using  the  same  physical  rationale  for  the  development  of  the  basic  flow  on  the 
elliptic  cone  one  may  attempt  to  neglect  the  lateral  and  wall-normal  velocity  components,  u  =  v  =  0. 

Contract  No.  F61775-00-WE069 


Inviscid  global  instability  of  flow  on  an  elliptic  cone 


13 


One  radical  simplification  from  a  numerical  point  of  view  is  offered  by  the  inviscid  limit  of  the 
equations  of  motion,  Re  — >  oo,  M  — >  0,  which  removes  both  these  sources  of  stiffness  in  the  system  to 
be  solved.  However,  this  resulting  system  is  a  coupled  problem  for  five  variables,  say  (p,  u,  v,  w.  T)t.  In 
this  respect,  the  issue  of  size  of  the  eigenvalue  problem  remains  unresolved.  Furthermore,  an  inviscid 
analysis  introduces  a  major  currently  unresolved  physical  issue,  that  of  the  extension  of  the  concept 
of  a  critical  layer  (introduced  and  unambiguously  defined  in  the  context  of  one-dimensional  instability 
analysis)  in  two  resolved  spatial  directions.  A  first  step  towards  resolving  the  issue  of  critical  layers  in 
two  spatial  dimensions  numerically  is  discussed  in  the  next  section. 

( b )  Oneinviscidlim.it 

The  only  means  by  which  an  inviscid  analysis  could  be  performed  at  a  lower  numerical  effort 
compared  with  the  various  forms  of  the  original  problem  is  the  case 

u  —  v  —  0,  w  /  0,  (2.27) 

which  results  in  the  compressible  generalisation  of  the  Rayleigh  equation  [25] 


Cp  + 


Vflx 

_  Px) 

2  f3wx 

Px  + 

7  Py 

_  Py) 

2  flWy 

Py  + 

p{flw-u)2 

At p 

p  ) 

(f3w  —  a>)_ 

At p 

P  ) 

(/ 3w  —  co)_ 

IP 

(2.28) 


where  C  =  d/dx2  +  d/dy2  —  ft2.  This  most  simple  form  in  which  the  compressible  two-dimensional 
eigenvalue  problem  may  be  recast  can  be  solved  as  a  cubic  equation  for  either  a  complex  temporal 
eigenvalue  to.  taking  (5  to  be  a  constant  real  wavenumber  parameter  in  the  z  spatial  direction,  or  for 
the  complex  spatial  eigenvalue  ft  taking  the  frequency  parameter  w  to  be  a  real  constant.  Notably 
there  is  no  overhead  when  solving  the  spatial  compared  with  the  temporal  two-dimensional  eigenvalue 
problem  (2.28). 

From  a  physical  point  of  view  the  implication  of  addressing  (2.28)  as  opposed  to  (2.19-2.23)  is 
that  the  three-dimensional  elliptic  cone  geometry  is  reduced  to  a  three-dimensional  elliptic  cylinder 
geometry  which  is  periodic  in  the  third  direction  with  a  periodicity  length  Lz  —  2ir/ ft  The  analysis  may 
be  performed  at  a  z  —  zq  —  const,  location  along  the  axis  of  rotation  of  the  cone  taking,  of  course,  as  a 
basic  flow  an  extract  from  the  correct  three-dimensional  basic  steady  state  already  calculated.  Having 
performed  the  analysis,  one  needs  to  compare  the  wavelength  of  any  globally  unstable  mode  predicted 
with  the  lengthscale  of  the  neglected  three-dimensionality  of  the  elliptic  cone  surface  downstream  of 
the  location  z  —  zq.  If  the  first  length  scale  is  much  smaller  than  the  second,  the  predicted  instability 
may  be  physically  realisable.  Otherwise,  either  or  both  of  the  approximations  (2.27)  and  an  inviscid 
analysis  must  be  abandoned  in  favour  of  viscous  analysis  based  on  numerical  solution  of  the  full  system 
(2.19-2.23).  At  this  point  in  time,  which  of  the  two  scenaria  will  prevail  is  unclear;  the  potential  savings 
realised  by  the  inviscid  approach  make  the  latter  the  obvious  candidate  to  commence  global  instability 
analyses  on  the  elliptic  cone.  In  the  next  chapter  we  develop  the  two-dimensional  grid  on  which  an 
inviscid  global  instability  analysis  may  be  performed  and  provide  the  elements  forming  the  basis  to 
derive  the  three-dimensional  extension  of  the  global  eigenvalue  problem  for  the  application  at  hand. 
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3.  On  frames  of  reference  on  the  elliptic  cone 


In  the  preceding  discussion  we  justified  interest  in  a  global  linear  analysis  based  on  the  two- 
dimensional  eigenvalue  problem  from  the  point  of  view  of  numerical  feasibility.  From  a  physical  point  of 
view  this  approach  can  be  justified  using  one  of  the  fundamental  assumptions  of  boundary-layer  theory, 
namely  that  the  dependence  of  the  basic  flow  quantities  on  the  downstream  direction  (to  be  defined 
shortly  in  the  context  of  the  elliptic  cone)  is  0(1/  Re)  weaker  than  that  on  the  wall-normal  direction, 
Re  denoting  the  free-stream  Reynolds  number  of  the  flow.  Experimentation  and  three-dimensional 
numerical  solutions  of  the  equations  of  motion  for  the  basic  flow  have  verified  this  assumption;  this 
property  is  also  to  be  found  in  the  low  angle-of-attack  solutions  of  the  Navier-Stokes  equations  pre¬ 
sented  in  [25]. 

Regarding  instability  analyses  information  on  three-dimensionality  comes  principally  from  exper¬ 
iments  and  direct  numerical  simulations  on  planar  and  axisymmetric  geometries  [17,  21].  One-  and 
quasi  two-dimensional  linear  and  nonlinear  instability  analyses  are  performed  at  a  particular  location 
on  the  surface  of  the  object  examined.  Homogeneity  of  space  in  the  case  of  planar  geometries  and 
axisymmetry  of  the  circular  cone  permit  introduction  of  a  harmonic  decomposition  of  disturbances 
in  the  spanwise  or  circumferential  direction,  respectively.  As  has  been  mentioned,  the  dependence  of 
the  basic  flow  on  the  streamwise  direction  is  neglected  and  a  harmonic  decomposition  is  introduced 
in  this  direction.  The  latter  assumption  is  made  also  in  the  context  of  a  global  analysis  based  on  the 
two-dimensional  eigenvalue  problem.  As  far  as  the  downstream  location  is  concerned,  the  arc-length 
defined  along  the  generators  and  measured  from  the  tip  of  the  elliptic  cone  is  the  natural  choice  of 
coordinate  along  which  the  dependence  of  the  basic  flow  can  be  neglected  in  comparison  with  that 
along  the  other  two  spatial  directions.  This  choice  is  straightforward  if  the  oncoming  flow  is  at  zero 
angle  of  attack;  at  increasingly  large  nonzero  values  of  the  angles  of  attack  and  yaw/bank  the  choice 
is  less  clear  since  the  three-dimensionality  of  the  flowfield  is  enhanced  by  nonzero  values  of  the  angle 
of  attack;  nevertheless,  we  adhere  to  the  arc-length  as  the  downstream  coordinate  in  these  cases  also. 
As  the  case  is  with  the  one-  and  quasi  two-dimensional  analyses  the  wall-normal  direction  defines  one 
essential  spatial  direction  on  which  information  is  required  for  a  global  linear  analysis.  On  the  other 
hand,  the  definition  of  the  third  spatial  direction  on  the  elliptic  cone  is  not  unique. 

(a)  The  elliptic  confocal  coordinate  system 

In  anticipation  of  solving  a  two-dimensional  global  linear  eigenvalue  problem,  in  which  on  grounds 
of  feasibility  the  dependence  of  the  flow  on  the  third  direction  is  neglected,  the  natural  choice  of 
coordinate  system  to  express  the  basic  flow  quantities  and  solve  the  eigenproblem  is  the  confocal 
elliptic  coordinate  system  O^rj,  schematically  depicted  in  figure  1.  A  conformal  transformation  in  the 
complex  plane,  discussed  in  detail  in  appendix  A,  relates  the  cartesian  coordinate  system  Oxy  to  O^rj 
through 


x  —  c  cosh  £  cos  y  (3.1) 

y  =  csinh£sin?7  (3.2) 

in  two  spatial  dimensions,  while  the  extension  to  three  dimensions  is  trivial  and  leads  to  the  well- 
known  elliptic  cylindrical  coordinate  system  The  inverse  transformation  is  given  by 


£  =  -  Rejcosh  1(x  +  iy)}, 

(3.3) 

r)  —  -  Im{cosh-1(a;  +  iy)}. 

(3.4) 
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The  coordinate  lines  of  this  orthogonal  system  are  confocal  ellipses  £  =  cnst.  and  hyperbolae 
77  —  cnst.  with  common  foci  at  x  —  ±c  —  ±\/a2  —  U2.  The  semiaxes  themselves  are 


a  =  ccosh£,  (3.5) 

&  =  csinh£,  (3.6) 

which  ensures  satisfaction  of  the  equation  of  the  ellipse  in  cartesian  coordinates.  The  semi-axes  a 
and  b  and  the  focal  distance  c  of  the  ellipse  are  used  to  define  fV!  —  cnst.,  the  ellipse  describing  the 
wall  of  the  elliptic  cone  at  a  particular  location  zq,  through 


.  1  a  +  b  a  +  b 

o  ln - T  =  ln - • 

2  a  —  b  c 


(3.7) 


Chain  rule  for  partial  differentiation  is  used  to  relate  derivatives  in  the  two  coordinate  systems, 


(  df/dt  f  df/dx  \ 

V  df/dv  )~J{  df/dy  J 

The  inverse  Jacobian  matrix  of  the  transformation  is 

j_ i  _  c  /  sinh£  cos 77  —  cosh£  sin 77 
~  h2  \  cosh£  sin  77  sinh£  cos  77 

where  h  is  the  metric  of  the  transformation 


\  -  —  1 

( j  1 

-h  ^ 

)  h2 ' 

v  h 

n  J 

h  —  —  hv  —  c\J  sinh2  f  + 


sin2  77. 


(3.8) 


(3.9) 


(3.10) 


The  gradient  and  Laplacian  operators,  appearing  in  the  two-dimensional  eigenvalue  problems,  have 
the  forms 


grad  /  =  {fe/hJv/h) 

and 

V2  /  =  fe/h2  +  fm/h2 


(3.11) 

(3.12) 


respectively. 

Desirable  properties  of  this  coordinate  system  are  its  strong  clustering  of  points  in  the  neighbourhood 
of  the  wall,  itself  being  a  coordinate  line,  such  that  no  numerically  introduced  discontinuities  are  to  be 
expected.  Indeed  it  has  been  known  for  some  time  [14]  that  when  using  high-order  numerical  methods 
for  the  approximation  of  a  differential  operator  in  the  neighbourhood  of  a  solid  boundary  care  has  to 
be  taken  that  the  boundary  itself  is  approximated  with  the  same  degree  of  accuracy.  Furthermore,  this 
natural  clustering  at  the  wall  is  expected  to  be  advantageous  in  resolving  the  boundary  layer  and  its 
instabilities.  In  the  far-field,  the  exponential  decay  of  the  density  of  the  coordinate  lines  f  —  cnst.  as 
£  — >  00  permits  taking  the  outflow  boundary  of  the  computational  domain  well  away  from  viscous 
regions  and  the  only  limitation  in  the  choice  of  a  value  for  ^  is  expected  to  be  the  shock  location. 

The  physical  implication  of  an  instability  analysis  based  on  the  elliptic  confocal  coordinate  system 
is  that  the  three-dimensionality  of  the  elliptic  cone  is  locally  neglected  in  the  following  sense.  The 
basic  flow  is  calculated  by  a  three-dimensional  steady  scheme  and  may  be  extracted  on  an  arbitrary 
two-dimensional  grid.  Using  the  elliptic  confocal  coordinate  system  one  selects  a  particular  location 
zq  along  the  cone  axis  and  extracts  the  three  velocity  components,  density  and  pressure  on  a  plane 
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normal  to  the  axis  Oz.  The  inverse  transformation  (3. 3-3. 4)  permits  describing  the  non-conservative 
variables  on  the  coordinate  system,  while  the  discussion  of  this  paragraph  forms  the  basis  for 
transforming  any  of  the  forms  of  the  two-dimensional  eigenvalue  problem  discussed  in  the  previous 
section  into  the  elliptic  confocal  coordinate  system  and  solving  it. 

In  so  doing,  one  approximates  the  elliptic  cone  surface  locally  by  an  elliptic  cylinder.  Potentially 
unstable  global  eigenmodes  to  be  delivered  by  the  analysis  will  be  characterised  by  a  periodicity 
length  Lz  —  2ir//3  corresponding  to  the  most  unstable  wavenumber  parameter  (i.  A  comparison  of 
the  wavelength  of  the  most  unstable  global  eigenmode  with  the  length  scale  locally  characterising  the 
departure  of  the  actual  elliptic  cone  from  the  assumed  elliptic  cylinder,  A z.  Essential  for  a  meaningful 
global  instability  analysis  from  a  physical  point  of  view  is  the  satisfaction  of  the  condition 

Lz  «  Az  (3.13) 

which  implies  that  the  neglected  three-dimensionality  of  the  elliptic  cone  does  not  affect  locally  the 
development  of  an  unstable  z— wise  periodic  global  eigenmode. 

On  account  of  the  underlying  analytic  conformal  transformation,  the  elliptic  confocal  coordinate 
system  is  the  simplest  and  most  elegant  means  of  performing  a  global  instability  analysis  on  the  elliptic 
cone.  However,  this  is  not  the  only  possible  approach  and  in  what  follows  we  discuss  briefly  the  building 
blocks  of  two  different  methodologies  towards  arriving  at  three-dimensional  global  instability  analysis 
equations.  In  following  this  approach  we  intend  to  provide  alternative  three-dimensional  coordinate 
systems,  the  coordinate  lines  of  which  are  by  construction  mutually  orthogonal.  This  is  the  case  for 
the  three-dimensional  extension  of  the  elliptic  confocal  coordinate  system,  namely  the  elliptic  cylinder, 
but  not  for  the  elliptic  cone.  The  following  two  paragraphs  sketch  two  alternative  approaches  to  arrive 
at  consistent  means  of  neglecting  the  third  spatial  direction,  while  retaining  the  three-dimensionality 
of  the  elliptic  cone  surface.  The  resulting  two-dimensional  surfaces  are  expected  to  be  different  than 
the  plane  defined  by  the  cone  base  on  which  the  elliptic  confocal  coordinate  system  has  been  discussed. 
As  a  consequence  performance  of  global  instability  analyses  on  such  three-dimensional  surfaces  would 
require  an  extension  of  the  approach  based  on  (2.1).  Such  an  extension  will  be  considered  in  future 
efforts.  Ultimately,  comparisons  of  the  global  instability  analysis  results  on  the  two  different  two- 
dimensional  planes  could  embed  the  known  physical  mechanisms  into  a  2  framework  and  shed  light 
at  potentially  existing  additional  physical  mechanisms  that  cannot  be  captured  by  the  simple  one¬ 
dimensional  linear  theory. 

( b )  A  body-fitted  three-dimensional  coordinate  system, 

Tensor  analysis  may  be  used  to  arrive  at  a  general  three-dimensional  body-fitted  system  to  describe 
the  equations  of  motion.  In  this  case  a  transformation  between  a  given  cartesian  frame  of  reference 
0(x,y,z)T  and  an  orthogonal  body-fitted  coordinate  system  0(£,ry,  #)T,  in  which  £  denotes  the  co¬ 
ordinate  from  the  elliptic  cone  tip  along  the  cone  surface,  rj  denotes  the  wall-normal  and  #  the  angle 
between  the  cone  axis  of  symmetry  and  the  line  0£,  can  be  constructed  by  simple  trigonometric 
relationships.  On  the  plane  x  —  xq  —  0,  schematically  presented  in  figure  2, 


z0  =  0,  (3.14) 

Vo  =  Vi  +  V2  =  (£  -  A£)  sin#  +  — 'q—  =  (£  -T) tan#)  sin#  +  — ,  (3.15) 

(3.16) 

where  the  elliptic  cone  semiangle  #  is  in  general  a  function  of  the  angle  <p.  shown  in  figure  3.  The 
equation  of  the  ellipse,  expressed  in  terms  of  its  eccentricity  e  =  c/a  and  the  parameter  n  —  l>2 /a.  a 
and  b  being  the  semiaxes  of  the  ellipse,  is 
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R(cf)  = 


K 

1  +  e  cos  <f> 


(3.17) 


In  three-dimensional  space,  £,  9  and  R  are  functions  of  the  angle  (f>.  as  shown  in  figure  4.  Simple 
trigonometrical  relations  permit  solving  for  the  independent  variable 


6  ((f))  —  arctan 


'  1  K 

_L  (1  +  e  cos  </)) 


Combining  (3.17)  and  (3.18)  it  follows  that 


(3.18) 


a<t>)  =  if2  + 


(1  +  e  cos  cf>) 


(3.19) 


such  that  the  body-fitted  coordinate  system  0(£,  q,  9)T  is  related  with  the  cartesian  system  Oxyz 
through 


x 

V 

z 


([< 

([« 


q  tan  9  sin  9  + 


q  tan  9 


sin  9  + 


q  tan  9)  cos  9. 


q 


cos  9 
9 

cos  9 


COS  <■/), 
sin  f/). 


(3.20) 

(3.21) 

(3.22) 


The  objective  now  becomes  to  express  the  three-dimensional  compressible  continuity,  Navier-Stokes 
and  energy  equations  in  the  body-fitted  coordinate  system.  This  is  currently  pursued,  however  at  lower 
priority  compared  with  the  current  main  focus,  which  is  to  obtain  two-dimensional  eigenvalue  problem 
results  on  the  two-dimensional  elliptic  confocal  coordinate  system  introduced  in  the  previous  section. 


(c)  Differential  geometry  concepts  for  the  derivation  of  the  instability  equations 

In  this  section  we  follow  an  alternative  little  known  methodology  in  order  to  arrive  at  the  desired 
two-dimensional  global  stability  equations  as  the  limit  of  those  expressed  in  a  three-dimensional  body- 
fitted  coordinate  system.  Ultimately,  the  present  approach  and  that  of  the  previous  section  should  yield 
identical  results.  In  view  of  the  non-trivial  nature  of  the  systems  of  equations  in  question,  following 
two  alternative  methodologies  is  useful  in  cross-verifying  the  results  of  the  respective  derivations. 

An  essential  property  of  the  geometry  in  question,  introduced  by  the  departure  of  the  elliptic  cone 
from  axisymmetry,  is  the  non-orthogonality  of  the  azimuthal  direction  and  the  plane  defined  by  the  arc- 
length  and  wall-normal  directions.  This  is  illustrated  by  the  introduction  of  a  body-fitted  coordinate 
system]'  defined  by  the  arc-length  coordinate  s,  the  wall-normal  coordinate  n  and  a  coordinate  9  normal 
to  the  plane  Osn.  This  system  is  shown  in  the  upper  part  of  figure  5;  also  shown  is  the  azimuthal 
coordinate  f>.  which  is  normal  to  the  Oyz  plane.  Along  the  (s,  n,  9)  directions  we  introduce  unit  vectors 
es,  en.  eg  and  e^.  The  Cartesian  frame  of  reference  (x.  y,  z)  on  which  the  basic  flow  calculations  have 
been  performed  in  [25]  and  the  newly  introduced  frames  of  reference  are  related  by 


x  —  cos  (f(sa  +  nbL ), 

(3.23) 

y  —  sin  cf)(sb  +  naL ), 

(3.24) 

z  —  sL  —  nab. 

(3.25) 

t  in  general  different  than  that  of  the  previous  section  -  different  notation  underlines  this  point 
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(3.27) 

(3.28) 

(3.29) 


In  the  limit  case  a  —  b  —  r,  which  corresponds  to  a  cone  with  circular  base  of  radius  r.  </>  =  0 
and  e<f)  —  eg.  However,  if  a  /  b  the  third  component  of  eg  is  nonzero  and  eg  defines  the  tangent 
vector  along  the  2-7r— periodic  line  L[b2  —  a2]  sin  <■/>  cos  e/>.  Also,  eg  is  by  construction  tangent  to  the  line 
s  —  cnst.  on  which  the  elliptic  cone  is  shown  to  terminate  in  the  lower  left  part  of  figure  5  and  to 
the  wall-normal  direction.  In  this  case  •  eg  ^  0,  a  fact  which  has  far-reaching  consequences  for  the 
global  instability  analysis  in  terms  of  the  definition  of  the  plane  on  which  the  global  analysis  will  be 
performed  and  hence  we  elaborate  on  the  coordinate  systems  defined  further. 

The  two  coordinate  systems  (es,  en,  e^)  and  (es,  en .  eg )  have  the  following  properties.  Both  reference 
systems  have  in  common  a  local  origin  which  is  placed  as  the  same  distance  from  the  cone  tip, 
a  property  which  may  be  interesting  in  the  subsequent  analysis  if  we  wish  to  describe  the  three- 
dimensionality  of  the  flow  in  terms  of  its  impact  on  local  instability  properties.  Specifically,  local 
linear  analyses  could  be  performed  by  monitoring  the  flow  at  equidistant  locations  from  the  tip  of 
the  elliptic  cone  and  resolving  the  wall-normal  direction  and  the  results  could  be  compared  with 
those  of  global  analyses  in  which  the  s  —  cnst.  line  defines  one  resolved  spatial  direction.  The  second 
common  property  of  the  two  coordinate  systems  is  that  they  encompass  one  direction  normal  to  the 
elliptic  cone  surface;  this  is  an  essential  characteristic  independently  of  the  type  of  global  or  local 
instability  analysis  to  be  performed.  The  discriminating  characteristic  of  the  two  coordinate  systems 
is  that  (es,en,e^)  is  non-orthogonal  by  contrast  to  (es,  en.  eg)  which  is  orthonormal.  While  the  last 
characteristic  is  appealing  in  the  derivation  of  the  global  instability  equations  on  the  OnO  plane,  use  of 
the  (es,en,e^,)  coordinate  system  could  be  advantageous  in  facilitating  comparisons  with  experiment 
in  which  a  light  sheet  is  used  to  visualise  instabilities  on  the  Oncf)  plane  [12]. 

A  third  coordinate  system  may  also  be  introduced  along  the  following  lines.  The  arc-length  direction 
Os  and  its  unit  vector  es  as  well  as  the  local  wall-normal  direction  On  and  its  unit  vector  en  are 
introduced  as  before.  The  third  coordinate  direction  is  defined  in  this  coordinate  system  numerically 
by  reference  to  the  basic  flow  data  of  Theofilis  [25] .  The  three-dimensional  flowfields  could  be  post- 
processed  to  recover  the  magnitude  and  direction  of  the  flow  perpendicular  to  the  Ons  plane,  which 
we  loosely  refer  to  as  the  crossflow  direction.  The  third  coordinate  could  then  be  defined  by  the 


Contract  No.  F61775-00-WE069 


Inviscid  global  instability  of  flow  on  an  elliptic  cone 


19 


requirement 

crossflow  magnitude  =  cnst.  (3.30) 

This  introduces  a  different  direction  Ox) ,  termed  the  crossflow  direction,  along  which  a  tangent 
vector  e,j  could  be  defined  numerically.  The  global  linear  analysis  can  the  be  performed  on  any  of  the 
three  coordinate  systems  (es,en,e^),  (es, en.eo)  or  (es, en. e,j)  by  resolving,  respectively,  the  planes 
One/),  On6  or  Onx).  The  advantage  of  performing  the  analysis  on  the  Onx)  plane  would  be  that  the 
effects  of  wall-normal  and  crossflow  direction  on  flow  instability  would  be  separated  in  the  the  global 
analysis. 

The  benefit  of  choosing  any  of  the  three  coordinate  systems  to  perform  the  global  instability  calcu¬ 
lations  could  be  assessed  by  reference  to  experiment.  However,  the  form  of  the  equations  of  motion  is 
obviously  affected  by  choosing  different  coordinate  systems  and  the  nontrivial  nature  of  the  derivations 
of  the  three-dimensional  compressible  equations  of  motion  in  arbitrary  coordinate  systems  warrants  a 
decision  at  this  stage  on  the  reference  system  to  be  used  in  the  subsequent  analysis.  As  an  example, 
we  use  Appendix  C  to  demonstrate  derivation  of  the  gradient  operator  using  a  little-known  methodol¬ 
ogy  of  differential  geometry  [11,  3],  which  delivers  transformation  of  equations  in  arbitrary  coordinate 
systems  including  those  defined  numerically,  such  as  the  (es,  en,  e,j )  frame  of  reference,  without  having 
to  resort  to  tensor  analysis.  The  result  in  the  (es,en,e$)  coordinate  system  is 


df 


]_d£ 

cot  ds 


+ 


LO 


<t> 


LOj 


cotco 


1  df 

Vde 


LO 


1  df 


(3.31) 


with  the  shorthand  notation  explained  in  Appendix  C.  The  differential  geometry  concepts  are  useful 
in  providing  alternative  ways  of  deriving  the  three-dimensional  compressible  Navier-Stokes  equations  in 
general  three-dimensional  frames  of  reference,  the  only  requirement  for  which  is  that  they  be  orthogonal 
at  each  point  in  space.  In  view  of  the  relative  simplicity  of  the  three-dimensional  shape  of  an  elliptic 
cone  surface,  which  will  form  the  basis  for  the  derivation  of  an  body-fitted  frame  of  reference,  we 
turn  to  deriving  the  equations  of  motion  using  classic  tensor  analysis  and  will  return  to  differential 
geometry  in  order  to  cross- validate  the  results  of  the  anticipated  extensive  calculations. 
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4.  The  numerical  solution  of  the  generalised  Rayleigh  equation 

In  view  of  the  transformation  (3. 1-3.2)  the  generalised  Rayleigh  equation  (2.28)  becomes  on  the 
elliptic  confocal  coordinate  system 


PH  +  Pr)ri 


h2/32p  + 
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fl  +  32 

'(el 

h2 
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fl  +  32 
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2  /3wv 
(t 3w  —  uj)_ 
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p  —  0. 


However,  j\  +  j\  —  h2 ,  such  that  one  finally  has  to  solve 


(4.1) 
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2  (3u>£ 
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p  —  0. 


(4.2) 


where  the  linear  operator  Ai  =  +  drm  —  h2/32.  In  a  manner  analogous  with  classic  one-dimensional 

linear  theory,  (4.2)  may  be  solved  either  iteratively  or  by  direct  means.  In  view  of  the  lack  of  any  prior 
physical  insight  into  global  linear  disturbances  in  the  application  at  hand,  a  direct  method  is  preferable 
on  account  of  the  access  to  the  full  eigenvalue  spectrum  that  it  provides. 

Either  the  temporal  or  the  spatial  form  of  the  eigenvalue  problem  may  be  solved  at  the  same  level 
of  numerical  effort  using  a  direct  method  since,  in  both  cases,  a  cubic  eigenvalue  problem  must  be 
solved.  In  its  temporal  form,  this  is 


TiPii  +  t2 Prm  +T3Pi  +  Tipv+T5p  =  uj  (T6%  +  Tjpm  +  T8p^  +  Tgp^  +  T10p) 

+  w2T„p 

+  w3T12p  (4.3) 

while  the  spatial  generalised  Rayleigh  equation  on  the  elliptic  confocal  coordinate  system  is 

Si  Pa  +  S2pnn  +  S3p£  +  SiPn  +  S5p  -  f3  (S6%  +  +  S8ps  +  S9p v  +  Swp) 

+  P2  Sup 

+  P3  Sup  (4.4) 

The  coefficients  appearing  in  (4. 3-4.4)  are  given  in  appendix  B. 

(a)  Spectral  collocation  for  the  inviscid  instability  problem 
Either  of  the  systems  (4.3)  or  (4.4)  may  be  discretised  by  spectral  collocation  means.  The  natural 
choices  for  basis  functions  are  a  member  of  the  Jacobi  family  of  polynomials;  here  Chebyshev  polyno¬ 
mials  have  been  chosen  to  discretise  the  f—  spatial  direction  [4],  The  natural  boundary  of  the  elliptic 
cone  surface  dictates  the  need  for  a  mapping  of  the  standard  Chebyshev  domain  Xj  —  cos  jn/N^  onto 
£  through 


—  2  xj  (£■ w  £oo)  +  2  (£■ w  £oo) :  j  —  0,  ■  ■  ■ ,  Ns 


(4.5) 
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where  is  given  by  (3.7)  and  £oo  represents  the  far-held  location  along  the  f  coordinate  where  the 
domain  is  truncated.  Chain  rule  is  used  to  account  for  the  metric  of  the  transformation  between  Xj  and 
analytically  and  arrive  at  the  expression  for  the  collocation  derivative  matrices,  V^.  for  the  first  and 
T>cc  for  the  second  derivative  in  the  direction  in  terms  of  their  standard  Chebyshev  counterparts 

cfandB?. 


vf  = 


®«  = 


{£w  £oo) 

4 

(iw  -  £oo)2 


Vr 


v'i. 


(4.6) 

(4.7) 


A  Fourier  expansion  is  used  to  discretise  the  2-7r— periodic  spatial  direction  rj.  Here  no  mapping  is 
necessary  and  one  may  use  the  standard  Fourier  points, 

Vj  —  2  j  7r  /Nv ,  j  =  (4.8) 

and  the  the  first  two  associated  Fourier  collocation  derivative  matrices  Vf  and  Vj  to  perform 
differentiation  along  rj. 


V^Vf,  (4.9) 

Vm  =  Vj.  (4.10) 

The  generalised  Rayleigh  equation  may  now  be  discretised  in  elliptic  confocal  coordinate  space  to 
arrive  at  its  temporal 


{Ti  V^  +  T2  V^  +  T3  V^  +  T4  Vv  +  T5  J}p  —  u>  {T@  V^  +  T7  V^  +  V^  +  Tg  Vv  +  T40  /} 

+  oj2  Tu  I P 

+  w3  T12 1  p  (4.11) 

or  spatial  form 

{Si  V^  +  S2  Vvv  +  S3  V^  +  1S4  Vv  +  S5 1}  p  —  fl  {Sq  V^  +  S7  Vvv  +  Ss  V^  +  Sg  Vv  +  5io  7} 

+  fl2SnIp 

+  ft3  S12I p  (4.12) 


which  may  be  readily  collocated. 

However,  before  either  equation  can  be  solved,  the  issue  of  critical  layer  in  two  spatial  dimensions 
arises,  which  must  be  resolved.  In  the  next  chapter  we  provide  one  solution  based  on  the  idea  of 
complex  collocation  calculation  grids.  This  is  extending  spectral  collocation  grids  of  the  type  (4.5) 
and  (4.8)  into  complex  space,  where  the  calculations  are  performed;  the  physical  quantities  of  interest 
are  extracted  as  the  real  part  of  the  results  obtained.  For  simplicity  of  the  demonstration  and  in  view 
of  the  abundance  of  comparison  data  existing  we  adhere  to  the  classic  one-dimensional  compressible 
linear  stability  eigenvalue  problem. 
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5.  On  complex  spectral  collocation  grids  for  inviscid  instability  analysis 

In  anticipation  of  an  inviscid  global  analysis  to  follow,  we  discuss  here  the  development  of  a  new 
methodology  based  on  spectral  collocation  methods  to  address  this  issue.  While  (real-grid)  spectral 
methods  have  long  been  applied  to  the  solution  of  the  basic  flow  problem  [20]  their  application  to  the 
resultant  inviscid  instability  analysis  has  proved  prohibitive  due  to  the  existence  of  the  critical  layers. 
Two  solutions  exist  in  order  to  regain  spectral  accuracy.  The  first  is  to  use  a  standard  Chebyshev 
Gauss-Lobatto  (real)  grid  [4]  but  Taylor-expand  the  basic  flow  around  the  critical  layer.  The  second, 
followed  herein,  is  to  extend  the  Zaat-Mack  technique  for  the  integration  of  the  instability  equations 
in  the  complex  plane  and  introduce  complex  collocation  grids,  forcing  a  spectral  method  to  integrate 
the  equations  on  a  complex  contour  suitably  adapted  to  avoid  the  critical  layer.  Use  of  this  second 
approach  was  first  suggested  by  Boyd  and  Christidis  [2]  and  further  investigated  by  Boyd  [1]  in  the 
context  of  atmospheric  and  hydrodynamic  instability  calculations.  These  works  were  extended  by  Gill 
and  Sneddon  [10],  who  gave  an  analytic  formula  for  optimizing  one  family  of  complex  (quadratic) 
maps. 

Here  we  validate  the  key  element  of  an  inviscid  global  linear  instability  analysis,  namely  complex 
spectral  collocation  calculation  grids,  by  introducing  a  novel  spectral  scheme  for  the  solution  of  the  one¬ 
dimensional  incompressible  inviscid  linear  eigenvalue  problem  in  planar  geometries  and  that  pertaining 
to  supersonic  flow  in  both  planar  and  curved  geometries.  The  first  problem  is  governed  by  the  classic 
Rayleigh  equation  while  compressibility  and  curvature  are  discussed  using  the  equations  of  Duck  [8] 
or  their  planar  limits.  We  employ  complex  spectral  collocation  grid  techniques  first  to  the  solution 
of  the  basic  flow  problem  and  subsequently  to  both  the  incompressible  and  the  compressible  inviscid 
linear  eigenvalue  problems.  The  accuracy  and  the  efficiency  of  the  proposed  algorithms  are  assessed  in 
comprehensive  comparisons  of  results  obtained  using  our  spectral  approach  against  the  work  of  Mack 
[17,  18]  and  the  finite-difference  algorithm  of  Duck  [8]. 

(a)  Linear  local  inviscid  instability  theoretical  background 

With  reference  to  the  conical  geometry  presented  by  [9],  introduction  of  the  local  theory  Ansatz  into 
the  full  Navier-Stokes,  energy  and  continuity  equations  results  in  a  sixth-order  system  for  the  linear 
disturbance  amplitude  functions;  further,  if  one  takes  the  0(e  -C  l)-terms  with  the  leading  order  in 
Reynolds  number  (i.e. ,  ignoring  viscous  terms),  the  sixth-order  system  can  be  reduced  to 


Vn  + 


C 


l  +  AC2  +  c/  Wo -13  ^MKWo-P) 


Wqv<p  _ 


1 p 


2a2 


1  + 


n% 


a2(l  +  X(2  +  (r/)2 


-Ml(W0-/3f 


ia2(Wo  —  /3)-^r  —  — 


7^’ 


(5.1) 

(5.2) 


in  the  independent  variable  r).  where  Wo  and  To  are  respectively  the  single  basic  flow  velocity  com¬ 
ponent  and  basic  flow  temperature  profiles,  A  =  0  corresponds  to  a  cylinder  and  A  /  0  generates 
families  of  conical  geometries,  (  is  the  scaled  downstream  coordinate,  <p  —  v\ /£  and  p  are  the  scaled 
disturbance  velocity  component  and  disturbance  pressure  amplitude  functions,  M0 0  is  the  free-stream 
Mach  number,  7  is  the  ratio  of  specific  heats,  n  ^  0  corresponds  to  non-axisymmetric  disturbances, 
a  —  a(.  co  —  to(  and  (3  —  c  (c  complex  and  a  real)  for  temporal  instability  calculations  while  [3  —  co/a 
(a  complex  and  co  real)  for  spatial  instability  calculations  and  i  =  y/—l.  The  boundary  conditions  are 
[9] 


<p  —  prj  —  0  on 
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V  ~  {Kn+l{v)  +  ^jn-l|W)}  as 


r]  — >■  oo, 


(5.4) 


~  ^oo^ioiTl1  ~  fl)Kn{v) 

p~  [i  -Mui-m1/2 


as 


r]  — >  oo, 


(5.5) 


where  f)  —  ±a[  1  —  M^(l  —  /3)2]1/2(1/^  +  AC  +  r/),  the  sign  of  which  is  chosen  such  that  the  real  part 
of  the  argument  is  positive  in  order  for  boundedness  as  y  — >  oo  to  be  ensured,  ip ^  is  a  constant  and 
ifri(zi)  denotes  the  modified  Bessel  function  of  order  n  and  argument  z\. 

To  achieve  the  planar  limits,  we  firstly  set  A  =  0  in  (5. 1-5.2)  and  subsequently  apply  the  limit  £  — >  0. 
corresponding  to  flow  in  planar  geometries.  This  yields 


Vr/ 


WoflP  _ 


1 P 


Wo -13  ^MKWo-fl) 
ia2(W0-fl) 

The  corresponding  boundary  conditions  are 


T0  -  A4l(W0  py 


To  7^' 
f  —  Pv  —  0  at  r]  —  0, 


and 


~  <Poo  exp[— Ctyjl  -  (1  -  T)\,. 

'np00'yM%0(l  -  0)  exp[— a\/l  -  (1  -  fl)2M‘^  rj\ 

p  ~ - ,  ,  = -  as  r;  — >  oo. 

Vl-(l~(3)2Ml 

Combining  equations  (5. 6-5. 7)  and  taking  the  limits  — >  0  and  To  — >  cast.,  results  in 

(w0-P)[r  ~2- 


tprjri  ®  P 


WoririV  —  0j 


(5.6) 

(5.7) 

(5.8) 

(5.9) 

(5.10) 

(5.11) 


which  is  the  classic  Rayleigh  equation  governing  inviscid  linear  instability  of  incompressible  boundary- 
layer  flow  in  planar  geometries. 


(b)  Numerical  Methods 

The  physical  range  of  the  type  of  boundary  layer  being  considered  in  what  follows  extends  from 
zero  at  the  solid  boundary  to  suitably  chosen  far-held  positions.  This  necessitates  the  use  of  mapping 
transformations  between  this  range  and  the  domain  upon  which  the  Chebyshev  spectral  collocation 
points  are  defined.  The  standard  (real)  Gauss-Lobatto  collocation  points 

Xj  =  cos^,  (j  =  0,...,TV),;  (5-12) 

form  the  basis  of  all  complex  grids  constructed  in  the  sequel.  The  first  of  the  complex  grids  employed 
(and  the  most  straightforward)  is  based  on  transforming  [—1, 1]  to  a  parabolic  contour  in  the  complex 
domain  [1],  There  are  two  ways  to  construct  a  complex  grid  which  passes  from  the  point  r)  —  0  to 
rj  —  Vmax-,  the  location  on  the  real  axis  where  the  calculation  domain  is  truncated. 

One  approach  is  to  apply  first  the  complex  quadratic  transformation  [1] 

Vj  -  xj  +  iSflx2  -  1),  i  =  a/— T,  (5.13) 

taking  the  parabola  that  cuts  the  real  axis  at  y  —  —  1  and  y  —  1.  Here  8\  is  a  parameter  whose  effect 
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on  the  resultant  solutions  has  been  discussed  by  Theofilis  et  al.  [13].  The  physical  domain  is  then 
mapped  onto  the  computational  domain  using  the  following  algebraic  transformation  [23] 


”> =  ‘ttttvs  (5'14) 

where  s  —  21  /r;niax  and  l  is  a  stretching  parameter.  The  grid  defined  by  (5.14)  is  parabolic  and  avoids 
the  critical  point/layer  by  passing  below  it. 

Conversely,  the  interval  [—1,1]  may  be  transformed  to  the  computational  interval  [0,  ymax\  (where 
Umax  —  Vmax )  using  the  algebraic  transformation 


Vj  =  1 


1~xj 

1  +  S  +  Xj  ’ 


(5.15) 


and  then  maps  it  onto  the  complex  domain,  passing  through  the  real  points  rj  —  0,ry  =  ijmax ,  using 


Vj  =  Vj  +  i  h 


(-1)ro(w  -<W) 

Vmax 


(5.16) 


where  m  is  an  integer  with  m  >2  and  82  is  a  parameter,  both  of  which  determine  the  characteristics  of 
the  integration  path;  as  m  increases  (assuming  that  82  <  0  and  remains  constant)  the  maximum  depth 
of  the  complex  integration  path  increases  and  its  minimum  shifts  towards  the  Re(rj)  —  0  axis,  thus 
enabling  the  handling  of  critical  points/layers  which  lie  close  to  rj  —  0.  This  is  useful  in  calculations 
of  inviscid  instability  of  incompressible  wall-bounded  flows.  On  the  other  hand,  by  changing  the  sign 
of  8 2,  one  can  integrate  above  or  below  the  real  axis,  for  82  >  0  or  82  <  0,  respectively.  Further, 
by  increasing  the  absolute  value  of  S 2  one  can  increase  the  maximum  depth/height  of  the  complex 
integration  path.  Note,  the  real  parts  of  the  complex  grid  points  (5.14)  and  (5.16)  constructed  using 
these  alternate  approaches  are  not  identical.  In  other  words,  the  real  grid  generated  by  setting  <$1  =  0 
is  not  the  same  as  the  one  obtained  for  82  —  0. 

The  second  complex  grid  considered  is  dependent  on  the  availability  of  an  estimate  of  the  critical 
layer  position.  If  such  an  estimate  exists,  then  the  grid  can  be  deformed  locally  into  the  complex  plane 
using  an  exponential  complex  mapping  such  as  that  proposed  by  Boyd  [1] 


yj  =  Xj  +  i 83  e-(*l-*o)3/^j  (5T7) 

where  8 2 .  cq .  xq  are  parameters.  This  mapping,  combined  with  the  algebraic  transformation  (5.14), 
allows  for  a  short  detour  around  the  critical  layer  while  hugging  the  real  axis  for  the  rest  of  the 
interval.  Alternatively,  the  algebraic  mapping  (5.15)  combined  with  the  transformation 


Vj  =  Vj  +  y\  (5.18) 

could  be  used.  Note,  this  complex  grid  may  result  in  serious  convergence  problems  in  the  Chebyshev 
polynomials  expansion  of  the  eigenfunctions  [1].  As  in  the  case  of  real  grid  calculations,  the  chain  rule  is 
used  to  derive  the  modified  collocation  derivative  matrices  [4]  which  are  employed  in  the  computations. 


(c)  Instability  results  in  planar  and  axisymmetric  geometries 
The  classic  Rayleigh  equation  is  considered  first.  This,  despite  its  apparent  simplicity,  provides  a 
good  test  for  the  proposed  algorithm  due  to  the  fact  that  for  incompressible  boundary-layer  flows  the 
critical  layer  is  typically  located  very  close  to  the  wall.  The  complex  grids  which  we  employ  need  to 
coincide  with  the  real  point  on  the  fixed  boundary,  where  the  boundary  conditions  are  applied.  Thus 
the  computational  grid  must  be  deviated  around  the  critical  layer  at  a  very  short  distance  away  from 
the  boundary.  In  what  follows  we  use  two  sets  of  complex  mappings  resulting  from  a  combination  of 
(5.15)  with  either  (5.16)  or  (5.18),  both  of  which  have  the  ability  to  divert  the  complex  grid  over  short 
distances. 
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The  first  test  case  considered  is  the  Blasius  boundary  layer  [17]  which  has  a  critical  layer  location 
close  to  the  solid  wall.  The  basic  flow  profile  Wq(t])  was  obtained  by  solving  the  Blasius  equation  on 
a  real  grid  or  directly  on  a  complex  grid.  If  the  former  case,  a  Taylor  expansion  was  used  to  calculate 
the  values  of  W ( rj )  on  the  complex  nodes.  The  Rayleigh  equation  was  discretised  on  the  complex 
grid  resulting  from  the  combination  of  (5.15)  and  (5.16)  or  (5.15)  and  (5.18),  using  the  appropriate 
collocation  derivative  matrices.  The  discrete  eigenvalue  problem  was  thus  formulated  as  a  generalized 
eigenvalue  problem  which  was  solved  using  the  QZ  algorithm.  Mack  [17]  computed  eigenvalues  with  an 
indented  integration  contour.  For  a  —  0.179  Mack  quotes  the  eigenvalue  u  —  0.05750554  —  i0.00657109 
while  the  present  spectral  method  yields  u  —  0.05750493  —  i0.00657192  when  the  exponential  mapping 
with  <5i  =  —3.0  is  employed,  and  u  —  0.05750493  —  10.0065719 1  when  the  polynomial  mapping  with 
§2  —  —  0.15,  to  =  5  is  used.  For  both  spectral  method  calculations  100  collocation  points  were  used 
with  Tjmax  —  ^  —  60. 

The  second  test  case  is  the  flat-plate  boundary  layer  in  compressible  flow.  It  is  well-known  that  as 
the  Mach  number  increases  the  critical  layer  moves  away  from  the  wall  outwards  towards  the  edge  of 
the  boundary  layer  [17].  Thus  for  compressible  calculations  it  is  easier  to  determine  a  suitable  complex 
grid  which  avoids  the  large  gradient  problems  in  the  critical  layer.  By  contrast  to  the  Blasius  boundary 
layer,  where  an  inviscid  linear  instability  analysis  is  of  largely  academic  interest,  linear  instability  of 
(viscous)  compressible  flow  on  both  flat-plate  and  curved  geometries  can  be  well  approximated  by  an 
inviscid  analysis.  Mode  I  and  II  inviscid  linear  instability  calculations  were  performed  for  supersonic 
boundary-layer  flow  over  a  flat  plate,  choosing  the  relevant  parameters  from  the  related  work  of  Duck 
[8]  who  monitored  the  planar  limit  of  his  equations  against  the  calculations  of  Mack  [18].  Basic  flow 
calculations  were  performed  on  the  complex  collocation  points  defined  by  the  combination  of  (5.15)  and 
(5.16)  with  m  —  5.  The  subsequent  instability  analysis  was  pursued  in  a  temporal  framework  in  order 
for  comparisons  with  the  aforementioned  works  to  be  possible.  The  analysis  is  based  on  equations  (5.6- 
5.7)  together  with  boundary  conditions  (5.8-5.10).  First,  we  considered  the  incompressible  planar  limit 
taking  M0 0  =  10-4  and  the  mapping  parameters  of  the  incompressible  calculations,  N  —  100,  rjmax  — 
l  =  50,  <5-2  =  1.5.  The  calculated  eigenvalue  at  a  —  0.179  is  co  —  ac  =  0.05750493  —  iO. 00657192, 
which  agrees  with  the  incompressible  result  in  all  decimal  places.  Next,  supersonic  inviscid  linear 
instability  calculations  were  performed  at  —  3.8  for  flow  over  an  insulated  flat-plate.  We  chose  for 
all  subsequent  calculations  ijmax  —  l  —  25  and  present  in  Table  1  the  combined  effect  of  grid  resolution 
and  variation  of  the  single  remaining  complex  grid  parameter  82  on  the  mode  II  eigenvalue  obtained 
at  a  —  0.365.  In  these  results  it  may  be  seen  that  a  real  grid  is  inadequate  to  cope  with  this  instability 
problem.  The  same  conclusion  may  be  drawn  for  calculations  performed  with  <52  >  0.15  at  modest 
resolutions,  while  calculations  at  <52  ~  0.05  deliver  the  most  accurate  results  at  low  resolutions.  Finally, 
the  results  of  Table  1  as  well  as  others  not  shown  here  show  that  grid  refinement  in  combination  with 
a  small  positive  value  of  <52,  corresponding  to  a  short  detour  into  the  complex  plane,  can  deliver  the 
converged  eigenvalue.  Taking  the  optimal  parameters  resulting  from  the  present  study,  N  —  64  and 
<52  =  0.05,  we  present  in  Figure  6  the  dependence  of  the  eigenvalue  c  on  the  wavenumber  parameter 
a  for  both  mode  I  and  mode  II  calculations.  Excellent  agreement  with  the  superimposed  relevant 
inviscid  instability  analysis  results  of  Duck[8]  may  be  seen.  At  maximum  growth  rate  conditions  the 
viscous  analysis  of  Mack  [18]  predicts  a  mode  II  instability  wave  that  is  about  10%  more  stable  than 
the  inviscid  result  of  Duck [8]  reproduced  herein  by  the  complex-grid  analysis. 

The  final  validation  calculation  addressed  solution  of  the  linear  instability  equations  for  supersonic 
flow  past  axisymmetric  bodies.  The  temporal  linear  instability  analysis  results  of  [9]  are  used  for 
comparisons  with  our  predictions.  Figure  7  displays  results  obtained  with  the  Runge-Kutta  technique 
of  Duck  and  Shaw  [9]  and  our  algorithm  for  the  modes  I  and  II  of  instability,  for  the  case  of  a  cone 
(A  =  1),  at  £  =  0.05,  Mqo  =  3.8  and  azimuthal  wavenumbers,  n,  as  shown.  The  complex  grid  used 
is  defined  as  a  combination  of  (5.13)  and  (5.14)  with  N  —  64,  <5i  =  —0.05 —l  —  25.  The 
graphical  agreement  of  the  results  is  very  satisfactory.  Next,  we  perform  spatial  calculations  for  the 
same  physical  parameters.  We  use  a  combination  of  (5.13)  and  (5.14)  and  an  iterative  algorithm  at 
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variable  resolutions  to  recover  the  leading  eigenvalue.  The  dependence  of  the  eigenvalues  determined 
by  the  spectral  collocation  scheme  on  the  complex  grid  mapping  parameters  is  shown  in  Table  2. 
The  bracket  of  <5i  values  in  which  we  were  able  to  obtain  accurate  results  is  rather  narrow;  erroneous 
results  were  obtained  when  the  mapping  parameter  exceeds  a  threshold  value.  For  this  set  of  physical 
parameters  the  threshold  was  found  to  be  8\  ~  —0.05.  However,  the  principal  result  of  this  section 
is  that  the  complex  spectral  collocation  technique  proposed  is  capable  of  recovering  inviscid  linear 
local  instability  analysis  results  and  thus  becomes  one  viable  candidate  to  address  the  corresponding 
two-dimensional  eigenvalue  problem.  The  results  of  this  effort  will  be  presented  in  due  course. 
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(a)  Validation  runs 

(i)  Poisson  equation  on  the  elliptic  confocal  coordinate  system 

Our  first  concern  has  been  with  the  demonstration  of  the  tools  employed  for  the  numerical  solution 
of  the  generalised  compressible  Rayleigh  equation  on  the  elliptic  confocal  coordinate  system.  At  the 
heart  of  this  problem  lies  the  ability  to  solve  the  appropriate  Poisson  equation 


%  +  <V,  +  h2p(fl:  rj)  /(£,??)  =  q(£,ri), 


(6.1) 


with  h  —  h(£,r])  defined  in  (3.10)  and  p(£,  ??),</(£,  ??)  being  infinitely  differentiable  functions  on  the 
elliptic  confocal  coordinate  system  O^ig.  Boundary  conditions  to  close  the  system,  which  are  relevant 
to  the  subsequent  instability  analysis,  are  homogeneous  Dirichlet  at  the  farfield  and  either  Dirichlet 
or  Neumann  at  the  wall  of  the  ellipse. 

Several  examples  have  been  devised  and  solved,  in  order  to  validate  ( i )  the  coupled  two-dimensional 
spectral  collocation  discretisation  used,  (ii)  the  imposed  boundary  conditions  and  (Hi)  the  associated 
problem  of  symmetries  expected  (as  opposed  to  imposed)  in  the  solution.  In  all  cases 


has  been  considered, 


q(€,v)  =  1 

while  three  cases  have  been  monitored  for  p{^rj), 


(6.2) 


Case  H  :  p{€,v)  —  C 

and 

sinh  £  sin2  2rj 

CaseS:  p(£,q)  =  ,  == 

cy  sinh2  £  +  sin2  rj 

resulting  in  symmetric  solutions  about  both  axes  of  the  ellipse,  and 


(6.3) 


(6.4) 


Case  A  : 


p(t,v)  = 


sinh  £  sin  2r] 


c\j  sinh2  £  +  sin2  r\ 


3’ 


(6.5) 


which  breaks  this  symmetry  in  the  solution,  while  it  retains  generalised  symmetry  about  the  centre 
of  the  integration  domain. 

Natural  periodic  boundary  conditions  are  imposed  by  the  Fourier  expansion  along  the  r]  coordinate 
direction,  while  along  the  £  coordinate  either  Dirichlet 


/(£  =  tw,ri)  =  0,  /(£  =  £oo,v)  =  0  (6.6) 

or  Neumann 

m  =  ^,v)  =  o,  m  =  z  oo,v)=o  (6.7) 

boundary  conditions  are  imposed,  fw  denoting  the  coordinate  line  (3.7)  which  defines  the  ellipse 
surface  and  foe  being  the  location  along  the  f  coordinate  where  the  integration  domain  is  truncated. 

(n)  Numerical  solutions 

A  single  configuration  of  an  aspect  ratio  4  elliptic  cone  with  unit  major  semiaxis  has  been  considered. 
A  cut  plane  parallel  to  the  cone  base  at  a  distance  zq  —  0.5  from  the  cone  base  results  in  an  ellipse 
with  geometric  parameters 
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a  =  0.5,  b  =  0.125,  £w  «  0.255,  (6.8) 

while  ^oo  =  1.5  has  been  taken  to  truncate  the  integration  domain.  Six  problems  are  then  solved, 
namely  H,  S  and  A,  employing  Dirichlet  or  Neumann  boundary  conditions  at  the  wall.  The  spectral 
collocation  discretisation  scheme  of  §  4  is  used  and  the  convergence  history  of  the  results  for  the  wall 
value  f(x  —  a.  y  —  0)  =  /(£  =  ,  rj  =  0),  obtained  from  numerical  solution  of  the  Neumann  problems, 

is  presented  in  table  3  to  within  six  significant  digits,  while  the  convergence  history  of  the  numerical 
solution  of  problems  S  and  A  is  presented  in  figure  10.  The  spatial  distribution  of  f(x,y)  for  both 
Dirichlet  and  Neumann  boundary  data  can  be  found  in  figures  11-13. 

A  key  conclusion  based  on  these  results  is  that  in  all  cases  studied  exponential  convergence  of  the 
solutions  of  (6.1)  has  been  obtained.  The  symmetries  of  the  solutions  are  recovered,  rather  than  being 
imposed,  by  the  numerical  approach  taken.  Furthermore,  when  the  function  sought  has  a  simple, 
essentially  one-dimensional  (in  £)  structure,  as  the  case  is  in  problem  H,  a  small  number  of  collocation 
points  suffices  to  solve  (6.1)  and  further  increases  in  resolution  are  unnecessary.  This  result  is  attributed 
to  the  convergence  properties  of  the  Fourier  expansion  in  rj  and  underlines  the  significance  of  the 
choice  of  the  natural  elliptic  confocal  coordinate  system  for  the  problem  at  hand.  The  low  resolution 
requirements  in  the  rj  direction  in  this  class  of  problems,  effectively  suggesting  that  the  available 
computing  power  can  be  almost  exclusively  devoted  to  resolution  of  the  £  spatial  direction,  is  significant 
in  terms  of  the  ability  of  the  two-dimensional  eigenvalue  problem  to  recover  results  of  classic  one¬ 
dimensional  linear  theory;  in  this  case  the  instability  mode  can  be  described  with  a  small  number  of 
Fourier  modes  accounting  for  resolution  of  the  geometry  in  the  lateral  (rj)  direction,  while  its  principal 
variation  is  along  the  wall-normal  £  direction. 

As  the  structure  of  the  sought  function  becomes  increasingly  more  complicated  higher  resolution 
is  necessary  compared  with  that  needed  for  problem  H,  although  in  both  problems  S  and  A  also 
exponential  convergence  has  been  obtained.  However,  to  the  degree  that  conclusions  on  the  spatial 
structure  of  the  sought  BiGlobal  instability  modes  on  the  elliptic  cone  can  be  drawn  upon  evidence 
provided  by  the  numerical  solutions  of  problems  A  and  S,  a  resolution  issue  for  the  inviscid  compressible 
eigenvalue  problem  may  arise  if  accuracy  of  the  eigenvalue  problem  results  beyond  three  to  four 
significant  digits  is  necessary,  for  instance  in  the  case  of  near-neutral  modes.  In  such  cases  a  combination 
of  the  numerical  methods  of  §  4  and  §  5  may  be  necessary. 

( b )  Instability  analysis  results  on  the  elliptic  confocal  coordinate  system 
(i)  Extraction  and  mapping  of  basic  flow  data  onto  the  O^y  plane 

The  numerical  solution  of  the  compressible  generalised  Rayleigh  equation  (2.28)  requires  basic  state 
results  on  a  plane  z  —  zq.  parallel  to  the  elliptic  cone  base,  as  schematically  depicted  on  figure  8. 
Simple  trigonometric  formulae  relate  the  semiaxes  a,  b  and  focal  distance  c  of  the  ellipse  on  the  cut 
plane  with  the  respective  quantities  at  the  base  of  the  elliptic  cone,  A,  B ,  C.  through 


a  b  c  zo 
A  ~  B  ~  C  ~  ~Z' 


(6.9) 


For  a  given  elliptic  cone  geometry  and  distance  zq  of  the  cut  plane  from  the  tip  of  the  cone,  the 
quantities  a,  b  and  c  can  be  calculated  in  this  manner.  In  turn,  with  known  semiaxes  of  the  ellipse  on 
the  cut  plane  (3. 1-3.2)  relate  the  cartesian  system  Oxyz  of  the  basic  flow  calculations  (Theofilis  [25]) 
with  the  elliptic  confocal  system  Of  tj  on  which  the  instability  analysis  is  to  be  performed. 

Basic  flow  data  are  transferred  from  the  steady  three-dimensional  unstructured-grid  results  using 
second-order  accurate  interpolation  onto  an  optimal  ’nearest-neighbour’  two-dimensional  tessellation  of 
the  Oxy  plane,  the  latter  generated  by  the  interpolation  software  which  supports  the  3D-visualisation 
package  used.  Since  the  basic  flow  results  have  been  obtained  using  a  second-order  accurate  finite- 
volume  method,  use  of  a  higher  order  interpolation  procedure  is  not  justified  and  one  has  to  rely  on 
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high  resolution  in  order  to  improve  the  quality  of  the  basic  state.  A  second  bivariate  interpolation 
procedure  generates  a  unique  Thiessen  triangulation  and  calculates  basic  flow  derivatives  on  the  new 
unstructured  grid  on  the  Oxy  plane.  The  result  is  then  used  to  perform  a  third  and  final,  also  second- 
order  accurate,  interpolation,  which  delivers  the  necessary  basic  flow  streamwise  velocity  component, 
density  and  pressure  on  the  elliptic  confocal  coordinate  system  Ofrj.  Derivatives  of  these  quantities 
with  respect  to  £  and  rj  are  then  calculated  using  the  collocation  derivative  matrices. 

Clearly,  this  sequence  of  interpolations  carries  the  potential  of  cumulative  error  being  introduced, 
so  alternative  approaches  pointing  to  future  research  directions  are  briefly  discussed  here.  Use  of  the 
same  (z— independent)  elliptic  confocal  coordinate  system  to  obtain  basic  flow  data  and  perform  the 
instability  analysis  would  eliminate  the  error  introduced  by  interpolations.  However,  the  implication 
would  be  that  flow  on  an  elliptic  cylinder,  rather  than  on  an  elliptic  cone,  would  be  monitored  by 
such  an  approach  and  information  on  flow  three-dimensionality,  currently  obtained  through  successive 
instability  analyses  at  different  zo  locations,  would  be  lost.  As  has  been  discussed  in  §  3,  a  three- 
dimensional  grid  constructed  of  elliptic  confocal  planes  normal  to  the  cone  axis  at  different  zq  locations 
is  not  an  orthonormal  frame  of  reference;  in  other  words,  on  account  of  the  three-dimensionality  of  the 
elliptic  cone  geometry  it  is  not  possible  to  construct  a  three-dimensional  orthonormal  frame  of  reference 
having  the  two-dimensional  plane  Oxy  on  which  the  present  analyses  are  performed  embedded  to  it. 

Returning  to  the  coordinate  system  of  figure  1,  a  related  issue  that  must  be  addressed  in  this  context 
is  the  extent  of  the  domain  on  which  basic  flow  data  are  interpolated/retained  along  the  £  direction. 
Although  it  is  in  principle  possible  to  define  the  cut  plane  for  the  instability  analysis  so  that  it  either 
incorporates  or  avoids  the  Mach  cone  (in  supersonic/hypersonic  flow),  this  was  not  done  in  anticipation 
of  potential  resolution  issues  in  the  instability  analysis  should  the  first  course  of  action  be  taken  or 
the  ambiguity  of  defining  an  outer  boundary  in  the  £  direction  in  the  region  between  the  cone  surface 
itself  and  the  shock  location.  Consequently,  the  shock  location  itself  was  chosen  as  the  boundary  at 
which  homogeneous  Dirichlet  boundary  conditions  were  imposed  on  the  disturbance  pressure  in  all 
subsequent  calculations.  At  the  wall  homogeneous  Neumann  conditions  were  imposed. 

(n)  Inviscid  eigendisturbances  on  the  elliptic  cone 

With  these  considerations  in  mind,  attention  has  been  focussed  on  three  configurations,  an  aspect 
ratio  3  elliptic  cone  in  subsonic  viscous  (M  —  0.5,  lie  —  103)  flow,  an  aspect  ratio  4  elliptic  cone  in 
supersonic  inviscid  (M  —  4)  and  an  aspect  ratio  2  elliptic  cone  in  hypersonic  viscous  (M  =  8,  Re  —  103) 
flow.  The  subsonic  calculations  have  been  performed  on  a  3D  unstructured  grid  encompassing  ~  2  x  106 
points  on  the  half  plane  x  >  0;  the  converged  (in  time)  steady  results  were  mirrored  with  respect  to 
the  plane  x  —  0,  while  in  the  second  and  third  cases  flow  over  the  entire  cone  has  been  solved.  The 
subsequent  instability  analysis  is  expected  to  deliver  an  answer,  in  the  form  of  sign  of  the  amplification 
rate,  as  to  whether  the  assumptions  of  symmetry  in  the  case  of  the  subsonic  flow  and  steadiness  in  both 
cases  studied  are  physically  realisable.  The  supersonic/hypersonic  calculations  have  been  performed  on 
unstructured  3D  grids  encompassing  >  2.5  x  106  points.  Respectively,  the  angles  of  attack  a  —  20°,  5° 
and  a  —  0°  have  been  considered,  while  in  all  cases  a  zero  yaw  angle  has  been  taken.  Monitoring  three 
different  basic  states  at  several  different  parameters  parameters  may  make  physical  interpretation  of 
the  respective  results  and  their  relation  to  one  another  a  rather  challenging  task.  However,  the  primary 
focus  of  this  work  has  been  demonstration  of  the  tools  necessary  for  the  analysis  from  a  numerical 
point  of  view;  systematic  investigations  of  inviscid  BiGlobal  instability  from  a  physical  point  of  view 
is  a  task  of  future  research. 

The  basic  states  utilised  for  the  analyses  at  the  extreme  Mach  number  conditions  are  shown  graph¬ 
ically  in  figures  14  and  15  [25].  In  line  with  the  assumptions  underlying  the  compressible  generalised 
Rayleigh  equation,  only  the  streamwise  basic  flow  velocity  component  w  was  retained  in  the  analyses 
(alongside  the  density  and  pressure  distributions).  This  assumption  is  increasingly  valid  as  the  Mach 
number  increases;  for  example,  in  the  supersonic  and  hypersonic  cases  solved,  the  peak  w  value  is  more 
than  one  order  of  magnitude  larger  than  that  of  v  and  several  orders  of  magnitude  larger  than  that  of 
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the  velocity  component  u  along  the  lateral  spatial  direction  x.  Correspondingly,  the  results  obtained 
are  expected  to  be  increasingly  relevant  to  the  three-dimensional  physical  problem  as  M  increases. 
In  all  cases  solved  the  spatial  concept  has  been  followed,  in  which  oo  is  taken  to  be  a  real  parameter 
denoting  frequency  and  the  eigenvalue  f3  —  6r  +  \(5\  is  sought;  in  a  manner  analogous  with  classic  linear 
theory,  fjT  —  Re{/3}  represents  a  wavenumber  along  the  homogeneous  z  spatial  direction  along  which 
the  disturbance  grows  or  decays  if  f3\  —  Im{/3}  <  0,  or  >  0,  respectively.  A  final  point  relevant  to  all 
calculations  performed  is  that  the  QZ  algorithm  has  been  used  for  the  recovery  of  the  full  eigenspec- 
trum.  While  this  algorithm  is  highly  inefficient  and  its  application  may  even  lead  to  concealment  of 
physically  relevant  modes  on  account  of  the  modest  attainable  resolution  in  comparison  with  Krylov 
subspace  iteration  methods  [26],  use  of  the  QZ  algorithm  was  deemed  necessary  at  this  stage,  since 
no  prior  guiding  information  on  potentially  interesting  areas  of  the  eigenvalue  spectrum  is  available. 
Consequently,  on  account  of  the  0((N^Nn):>’)  CPU  time  scaling  of  the  QZ  algorithm,  only  modestly 
high  resolutions,  up  to  x  Nv  —  40  x  40  have  been  used. 

Figure  16  presents  the  spatial  eigenspectrum  obtained  at  M  —  0.5,  to  —  10-2.  The  calculations  are 
are  inconsistent  for  most  travelling  disturbances  and  only  stationary  neutral  pressure  eigenmodes  could 
be  converged  at  these  conditions  using  up  to  30  x  30  collocation  points;  higher  resolution  increased 
the  number  of  converged  stationary  disturbances  but  did  not  yield  converged  travelling  mode  results. 
Three  stationary  neutral  eigenmodes  are  shown  in  figures  17-19.  The  most  striking  feature  is  the 
symmetry  about  x  —  0  in  all  three  results,  alongside  the  absence  of  symmetry  in  the  azimuthal  £ 
direction  in  the  first  two  eigenmodes.  The  symmetry  about  x  —  0  is  a  result  of  that  of  the  basic  state, 
while  the  absence  of  azimuthal  symmetry  implies  that  the  respective  eigenmodes  cannot  be  obtained 
by  application  of  the  simplified  linear  analysis  concept  in  which  the  Ansatz 


q  =  q(£,f/)  +eq(0  exp  [*(A»y  +  pz  -  tot) 


(6.10) 


replaces  (2.1).  The  present  results  are  the  first  indication  that  employing  (6.10)  inhibits  large- 
scale  non-azimuthally-periodic  perturbations,  such  as  those  shown  in  figures  17-18,  from  revealing 
themselves  in  the  eigenvalue  spectrum.  For  such  non-azimuthally-periodic  disturbances  (2.1)  must 
be  used.  Focussing  on  the  result  of  figure  17  one  notes  that  this  eigenmode  possesses  a  structure 
which  is  resolved  to  (qualitatively)  the  same  good  degree  as  those  of  the  model  Poisson  problems 
shown  in  the  earlier  section.  The  pressure  eigenfunction  and  its  derivatives  are  smoothly  reduced 
to  zero  at  the  outer  boundary,  while  most  of  the  activity  takes  place  in  the  immediate  vicinity  of 
the  cone  surface.  The  differential  pressure  on  the  windward  and  leeward  sides  of  the  cone  is  clearly 
visible  in  the  results.  Further,  the  77— derivative  of  this  result  bears  a  striking  similarity  with  the 
spatial  structure  of  the  basic  velocity  profile.  Qualitatively  analogous  conclusions  may  be  drawn  from 
the  results  of  figure  18,  although  in  this  case  the  locations  of  primary  vortex  separation  is  clearly 
visible  in  the  results,  as  stagnation  points  in  the  disturbance  pressure  signal.  This  result  points  to 
the  potential  of  utilising  BiGlobal  instability  analysis  to  provide  an  alternative  framework  for  the 
well-studies  vortex  breakdown  problem  on  delta  wings,  a  problem  which  has  received  a  large  amount 
of  attention  using  (6.10)  but  virtually  none  by  application  of  the  more  appropriate,  in  our  opinion, 
decomposition  (2.1).  Such  a  study  could  be  pursues  by  systematically  departing  from  the  classic 
circular  cone  geometry  [21]  to  progress  towards  a  delta  wing  configuration  via  a  sequence  of  elliptic 
cones  at  aspect  ratios  AR  G  (l,oc)  .  Finally,  the  result  of  figure  19  is  representative  of  a  class  of 
pressure  eigenmodes  also  found  in  the  simulations,  which  may  be  considered  as  azimuthally-periodic 
and,  as  such,  might  be  able  to  be  addressed  by  the  simplified  concept  (6.10).  However,  although  using 
(6.10)  at  a  single  (£w,  rj)— location  may  be  orders-of-magnitude  less  CPU-time-consuming  than  use 
of  BiGlobal  instability  analysis,  the  simplified  linear  theory  concept  must  be  repeated  at  a  sufficient 
number  of  rj— locations  along  the  wall  £w  in  order  for  a  picture  of  the  instability  mode  along  the  entire 
rj— coordinate  to  be  formed,  such  that  BiGlobal  analysis  may  be  competitive  even  for  this  subclass  of 
disturbances.  All  conclusions  put  forward  in  the  subsonic  case  must  be  viewed  in  the  framework  of 
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the  certainly  questionable  in  subsonic  flow  simplifying  assumption  u  —  v  —  0.  The  results  obtained 
may  be  viewed  as  approximations  of  the  BiGlobal  eigenmodes  which  are  expected  to  be  obtained 
if  the  full  viscous  problem  is  solved.  Nevertheless,  solution  of  the  latter  problem  is  not  expected  to 
discredit  the  main  result  obtained  herein,  namely  the  existence  of  large-scale  non-azimuthally  periodic 
eigendisturbances.  As  mentioned  before,  the  assumptions  of  the  compressible  generalised  Rayleigh 
equation  are  increasingly  valid  as  the  flow  Mach  number  increases,  and  we  next  turn  to  the  supersonic 
and  hypersonic  cases  addressed. 

Figure  20  presents  the  spatial  eigenspectrum  obtained  at  M  —  4.0,  lo  —  10-2.  By  contrast  to  the 
spectrum  of  figure  16  here  both  stationary  and  travelling  disturbances  could  be  converged  using  up  to 
30  x  30  collocation  points.  Attention  is  focussed  here  on  four  families  of  eigenmodes,  indicated  as  A, 
B,  C  and  D  in  figure  20,  which  at  least  qualitatively  have  converged  and  point  to  different  classes  of 
eigendisturbances;  representative  members  of  these  families  are  respectively  presented  in  figures  21-  24. 
In  comparison  with  the  respective  subsonic  modes,  all  supersonic  pressure  eigendisturbances  recovered 
were  found  to  be  more  compact,  extending  over  a  narrower  range  that  their  subsonic  counterparts, 
in  line  with  analogous  experience  with  modes  of  classic  linear  theory.  Members  of  family  A,  one  of 
which  is  shown  in  figure  21,  extend  only  on  the  windward  face  of  the  elliptic  cone  and  may  have  a 
single-  (as  that  shown  in  fig.  21)  or  multiple- mode  structure  along  the  azimuthal  direction.  Modes 
pertaining  to  family  B,  on  the  other  hand,  one  representative  member  of  which  is  shown  in  fig.  22, 
are  the  conceptual  leeward  analogon  of  the  windward  modes  A  and  exhibit  are  in  all  other  respects 
similar  characteristics  with  the  A  modes,  namely  a  compact  structure  in  the  neighbourhood  of  the 
body  surface  that  quickly  decays  to  zero  as  the  free-stream  is  approached.  Modes  pertaining  to  either 
of  families  A  or  B  might  be  well  approximated  by  classic  analysis  based  on  (6.10),  provided  the 
latter  is  employed  only  on  the  appropriate  face  of  the  cone.  Neutral  modes,  members  of  families  C 
and  D,  possess  a  structure  extending  over  the  entire  integration  domain;  two  examples  are  shown  in 
figures  23  and  24.  A  conclusion  with  far-reaching  implications  regarding  results  of  family  C  is  that 
BiGlobal  instability  analysis  may  have  a  useful  role  to  play  in  identifying  travelling  neutral  pressure 
eigenmodes  which  play  a  central  role  in  the  field  of  computational  aeroacoustics  (CAA).  Further 
research  to  quantify  this  conjecture,  which  could  form  one  potential  extension  of  the  present  work, 
is  required  in  this  area;  if  this  conjecture  is  confirmed,  orders-of-magnitude  computational  savings 
may  be  realised  in  CAA  by  applying  BiGlobal  instability  analysis  instead  of  the  currently  almost 
exclusively  used  DNS  methodology  [5].  Before  making  final  statements  in  this  respect,  the  influence  of 
the  domain  extent  on  the  instability  characteristics  of  the  neutral  travelling  pressure  eigenmodes  must 
be  assessed,  for  example  by  comparison  with  DNS  studies  of  the  same  geometry  [29,  30].  One  question 
that  may  be  posed  in  this  respect  is  whether  the  structure  of  family  C  BiGlobal  eigenmodes  is  a  result 
of  the  imposed  homogeneous  Dirichlet  boundary  conditions  at  the  shock  location,  and  consequent 
reflections  within  the  resolved  domain,  or  whether  imposition  of  (in)homogeneous  Neumann  far-held 
pressure  boundary  conditions  would  preserve/alter  the  recovered  structure,  for  example  in  terms  of 
the  clearly  identifiable  wavelengths.  Finally,  in  line  with  the  analogous  result  obtained  in  the  subsonic 
case,  stationary  neutral  eigenmodes  are  also  found  in  the  BiGlobal  spectrum,  here  identified  as  family 
D,  a  typical  member  of  which  is  shown  in  figure  24.  Both  the  different  levels  of  disturbance  pressure  on 
the  wind-  and  lee  side  of  the  how  as  well  as  the  slightly  non-symmetric  nature  of  this  mode  can  be  seen. 
The  latter  may  be  attributed  either  to  poor  spatial  resolution  of  the  basic  state,  although  it  is  likely 
that  the  symmetric  about  x  —  0  counterpart  of  this  mode  may  also  be  found  in  the  eigenspectrum. 
Another  issue  relevant  to  the  existence  of  inviscid  neutral  modes,  either  stationary  or  travelling,  has 
been  introduced  in  §  5,  namely  the  ambiguity  in  the  dehnition  of  a  critical  ’layer’  concept  in  the 
solution  of  the  inviscid  two-dimensional  eigenvalue  problem.  Even  if  one  is  satished  from  a  physical 
point  of  view  that  the  present  results  are  plausible  solutions  of  the  BiGlobal  eigenvalue  problem,  the 
neutral  results  obtained  herein  should  be  contrasted  against  of  solutions  of  (2.28)  employing  complex 
tensor-product  grids  before  proclaimed  as  being  physically  relevant,  which  defines  another  potentially 
interesting  extension  of  the  present  work. 
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Turning  to  the  highest  Mach  number  value  monitored,  one  notes  that  on  account  of  a  zero  angle  of 
attack  the  basic  state  is  symmetric  about  both  x  —  0  and  y  —  0,  which  results  in  analogous  symmetries 
being  preserved  in  the  BiGlobal  eigenvalue  problem  as  well;  a  typical  result  is  shown  in  figure  25.  The 
far-held  location  ^  along  the  coordinate  where  homogeneous  Dirichlet  boundary  conditions  on 
the  disturbance  pressure  have  been  imposed,  namely  ^  <  £s,  £,s(rl)  denoting  the  shock  location, 
has  been  kept  the  same  as  in  the  M  —  4  calculations.  The  main  qualitative  analogy  between  the 
BiGlobal  eigenspectra  at  M  —  4  and  M  —  8  is  that  (alongside  stationary  and  neutral  modes)  families 
of  travelling  modes  having  analogous  structure  to  families  A  and  B  are  identified  at  M  —  8  too.  To 
the  extent  that  comparisons  are  permitted  without  having  established  that  the  modes  in  question  at 
M  —  8  can  be  continuously  obtained  by  systematically  increasing  the  how  Mach  number,  the  travelling 
modes  of  a  specihc  family  at  M  —  8  are  found  to  be  more  stable  than  their  counterparts  at  M  —  4,  a 
result  which  is  in  line  with  the  predictions  of  classic  linear  theory  based  on  (6.10)  on  wave- like  linear 
disturbances. 

In  view  of  the  qualitatively  analogous  results  obtained  in  the  M  —  8  spectrum  compared  with  the 
M  —  4  case,  no  further  discussion  of  the  BiGlobal  eigenspectra  is  included  here.  Instead,  a  different 
path  in  the  multiparametric  space  at  hand  has  been  followed,  namely  investigation  of  the  effect  of 
extending  the  domain  on  which  the  BiGlobal  eigenvalue  problem  has  been  solved  well  beyond  the 
shock  location.  Motivation  for  this  undertaking  is  provided  by  the  different  significance  of  the  concepts 
of  ’wall-normal  distance’  and  ’far-held’  when  employing  (6.10)  or  (2.1).  In  the  case  of  classic  linear 
theory  a  characteristic  length  scale  in  experimentation  and  calculations  is  provided  by  the  thickness 
of  the  boundary  layer  developing  on  the  elliptic  cone.  By  contrast,  the  BiGlobal  instability  analysis 
results  obtained  so  far  at  the  lower  Mach  numbers  point  to  the  existence  of  large-scale  disturbances 
which  scale  with  the  elliptic  cone  semiaxes.  In  this  framework  the  shock  location  is  too  close  to  the  cone 
surface  to  be  considered  as  an  appropriate  ’far-held’  position  at  which  homogeneous  Dirichlet  boundary 
conditions  on  the  disturbance  pressure  may  be  imposed.  On  the  other  hand,  enlarging  the  extent  of  the 
integration  domain  approximately  by  a  factor  four  along  the  £  coordinate,  in  combination  with  the  use 
of  an  unmapped  grid  in  the  direction  which  relies  on  the  natural  clustering  of  the  elliptic  confocal 
system  to  resolve  near-wall  structures,  may  result  in  loss  of  resolution  adjacent  to  the  cone  surface 
and,  possibly,  inability  to  resolve  wave-like  disturbances  scaling  with  the  boundary  layer  thickness. 
It  is  precisely  the  objective  of  this  part  of  the  present  investigation  to  determine  whether  certain 
instabilities  exist  that  are  insensitive  to  this  domain  change  and  determine  their  spatial  structure. 

The  most  interesting  results  of  the  latter  investigation  have  been  the  following.  First,  analogously 
with  the  narrower  domain  computations,  stationary  and  travelling  neutral  modes  exist,  while  a  strong 
modification  of  the  location,  in  (/3r,  jH\)  parameter  space,  of  the  travelling  eigendisturbances  may  be 
observed  on  account  of  the  wider  integration  domain.  The  spectra  at  the  resolutions  utilised  which, 
despite  not  having  converged  yet,  demonstrate  qualitatively  consistent  structure  of  their  branches,  are 
shown  in  figure  26.  Second,  the  spatial  structure  of  stationary  and  travelling  neutral  modes,  which 
are  marginally  resolved  at  the  resolutions  utilised,  points  to  the  fact  that  these  modes  are  sensitive  to 
the  location  of  the  imposition  of  the  boundary  condition.  A  typical  result  may  be  found  in  figure  27, 
where  it  can  be  seen  that  no  physical  justification  exists  as  to  why  homogeneous  Dirichlet  boundary 
conditions  should  be  imposed  on  the  disturbance  pressure  at  the  chosen  value;  a  different  type  of 
boundary  condition  is  necessary  to  capture  these  modes  correctly.  Third,  by  contrast,  homogeneous 
Dirichlet  pressure  boundary  conditions  appear  to  be  adequate  for  travelling  disturbances;  members  of 
two  different  classes  (i.e.  branches  of  travelling  modes  in  the  eigenspectrum  of  figure  26)  are  shown 
in  figures  28  and  29.  Interesting  in  these  results  are,  at  least,  two  facts.  On  the  one  hand,  the  shock 
location  (which  may  be  inferred  by  comparison  with  the  result  of  figure  25)  does  not  appear  to  interfere 
with  /  influence  the  spatial  structure  of  such  eigendisturbances  which  are  smoothly  reduced  to  zero 
in  the  far-field.  On  the  other  hand,  these  two  classes  of  modes  are  essentially  different  in  one  aspect, 
namely  their  rj— dependence:  While  the  mode  shown  in  figure  28  could  be  expected  to  be  recovered  by 
application  of  (6.10)  at  appropriately  higher  resolution,  that  shown  in  figure  29  appears  to  fall  outside 
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the  scope  of  the  simplified  linear  analysis  and  require  the  present  BiGlobal  concept  for  its  study.  In 
other  words,  within  the  framework  of  the  assumptions  of  the  present  approach,  there  appear  to  exist 
in  the  elliptic  cone  linearly  amplified  eigendisturbances  which  are  localised  in  the  neighbourhood  of 
the  cone  itself,  scale  with  an  0(1)  geometrical  length  of  the  body  in  question  and  are  beyond  the 
scope  of  the  simplified  linear  theory.  Further  study  to  quantify  this  discovery,  which  is  consistent  at 
all  Mach  numbers  examined  and  could  utilise  existing  experimental  data  or  be  performed  in  parallel 
with  new  experimental  efforts,  is  essential  for  progress  in  the  theoretical  understanding  of  instability 
mechanisms  operative  in  this  technologically  key  flow  configuration. 
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7.  Discussion 


The  present  effort  has  commenced  with  the  global  instability  analyses  of  the  steady-state  basic  flows 
on  elliptic  cones  obtained  in  [25].  The  three-dimensionality  of  the  elliptic  cone  introduces  a  freedom 
in  the  choice  of  coordinate  systems  which  may  be  utilised  to  isolate  different  factors  affecting  flow 
instability  in  a  natural  manner  by  appropriate  choice  of  the  coordinate  system  on  which  the  analysis 
is  performed.  Three  different  alternatives  have  been  introduced,  having  distinct  properties  the  merits 
of  which  have  been  assessed.  On  grounds  of  numerical  efficiency  it  has  been  proposed  to  conduct  the 
analysis  on  the  elliptic  confocal  coordinate  system,  which  is  the  two-dimensional  analogon  of  the  elliptic 
cylindrical  system  in  three-dimensional  space.  This  mapping  has  (at  least)  three  advantages;  first,  it 
takes  into  account  the  essential  elliptical  nature  of  the  geometry  surface  in  a  natural  manner.  Second, 
it  naturally  clusters  points  in  the  boundary  layer  while  the  grid  density  decreases  exponentially  in 
the  far-field,  thus  permitting  placing  the  outflow  boundary  at  a  large  distance  from  the  body  surface 
without  wasting  gridpoints  in  that  region.  Third,  the  analytical  nature  of  the  mapping  permits  the 
calculation  of  metrics  of  the  transformation  without  introduction  of  numerical  error.  A  potential 
disadvantage  of  the  elliptic  confocal  coordinate  system  is  that  three-dimensionality  of  the  flow  in  the 
streamwise  direction  is  neglected;  however,  this  is  consistent  with  the  assumption  of  two-dimensional 
global  linear  theory  which  neglects  the  dependence  of  the  flow  on  the  third  spatial  direction  along  the 
elliptic  cone  axis,  i.e.  derivatives  of  flow  quantities  along  this  direction  are  neglected.  In  this  manner 
no  additional  physical  information  is  expected  to  be  lost  on  account  of  the  choice  of  the  coordinate 
system  compared  with  that  implied  by  two-dimensional  BiGlobal  analysis. 

The  two-dimensional  viscous  compressible  partial-derivative  eigenvalue  problem  has  been  derived 
in  a  cartesian  coordinate  system.  Potential  simplifications  have  been  discussed  one  of  which,  the  gen¬ 
eralised  Rayleigh  equation,  has  been  expressed  in  the  elliptic  confocal  coordinate  system  as  either  a 
temporal  or  a  spatial  eigenvalue  problem;  numerical  prescriptions  have  been  provided  for  its  colloca¬ 
tion.  However,  if  a  global  inviscid  instability  analysis  methodology  is  followed,  the  common  approach 
of  shooting/iteration  used  in  local  inviscid  analysis  is  inapplicable  and  the  analysis  must  be  performed 
by  solving  the  partial  derivative  eigenvalue  problem  as  a  nonsymmetric  generalised  matrix  eigenvalue 
problem.  This  may  lead  to  inaccurate  results  on  account  of  poor  resolution  of  the  critical  layer(s)  if 
computations  are  performed  on  a  real  grid  and  weakly-amplified/damped  or  (near-)  neutral  modes 
are  sought.  In  order  to  circumvent  this  problem  an  algorithm  has  been  presented  which  poses  the 
inviscid  local  eigenvalue  problem  as  a  matrix  eigenvalue  problem  defined  on  complex  spectral  colloca¬ 
tion  calculation  grids.  Since  this  algorithm  could  become  the  cornerstone  of  a  potential  global  inviscid 
analysis,  its  accuracy  and  robustness  has  been  demonstrated  in  several  validation  calculations  of  the 
local  inviscid  eigenvalue  problem  in  planar  and  curved  geometries. 

The  numerical  tools  proposed  for  the  solution  of  the  compressible  generalised  Rayleigh  equation  have 
been  validated  on  model  Poisson  problems,  devised  to  conform  with  the  boundary  conditions  and  mimic 
the  spatial  structure  of  the  expected  eigendisturbances;  exponential  convergence  of  the  model  problem 
solutions  has  been  demonstrated.  Attention  has  then  been  focussed  on  numerical  aspects  of  the  analysis 
of  three  elliptic  cone  configurations  at  a  variery  of  parameters/conditions  in  order  for  experience  with 
the  different  facets  of  this  novel  problem  to  be  accumulated.  Results  obtained  at  all  Reynolds  and  Mach 
numbers  studied  point  to  the  existence  of  large-scale  instabilities,  which  scale  with  a  geometric  length 
of  the  respective  elliptic  cone  and  not  with  the  thickness  of  the  boundary  layer  developing  on  its  surface; 
indeed,  the  persistence  of  large-scale  eigendisturbances  superimposed  upon  inviscid  basic  states  a-priori 
excludes  the  possibility  that  such  disturbances  can  be  associated  with  Tollmien-Schlichting  instability 
on  the  cone  surface.  Nevertheless,  numerical  evidence  (in  the  form  of  well-resolved  azimuthally-periodic 
eigenmodes)  has  been  obtained  that  the  latter  type  of  disturbance  may  also  be  resolvable  if  the 
integration  domain  is  restricted  within  the  Mach  cone  (in  supersonic/hypersonic  flow).  A  persistent 
finding  of  all  analyses  has  been  the  appearance  in  the  spectrum  of  neutral  stationary  and  travelling 
disturbances.  Although  several  such  eigenmodes  are  converged,  the  complex  grid  technique  presented 
must  be  used  before  neutral  or  marginally  amplified/damped  modes  are  proclaimed  to  be  physically 
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relevant.  On  the  other  hand,  travelling  disturbances  with  0(1)  amplification  rates  have  been  identified, 
which  can  be  reasonably  well  resolved  on  real  grids:  while  the  maximally  attained  resolution  at  which 
the  model  Poisson  problems  have  been  converged  to  within  several  significant  places  cannot  be  achieved 
on  present-day  hardware  when  using  the  QZ  algorithm  for  the  eigenvalue  problem,  specific  parts  of 
the  eigenspectrum  could  be  monitored  with  high  resolution  when  using  Krylov  subspace  iteration  (e.g. 
[26]). 

Several  questions  are  left  open  and  new  issues  have  been  generated  by  the  present  effort.  The 
apparent  ability  of  the  methodology  to  address  neutrally-stable  travelling  pressure  eigendisturbances 
points  to  a  potential  utilisation  of  BiGlobal  instability  analysis  in  a  CAA  context;  such  efforts  should 
be  accompanied  by  theoretical  considerations  regarding  the  concept  of  critical  ’layer’  of  the  two- 
dimensional  eigenvalue  problem  and  the  complex-grid  algorithm  presented  herein.  Subsonic  flow  has 
been  found  easiest  to  resolve  numerically,  albeit  that  the  restrictive  assumptions  of  the  generalised 
Rayleigh  equation  may  not  permit  studying  such  flows  from  a  physical  point  of  view  when  using 
this  type  of  inviscid  analysis.  Nevertheless,  the  reasonably  accurate  description  of  the  area  of  flow 
separation  from  an  elliptic  cone  surface  at  an  angle  of  attack  points  to  the  potential  of  BiGlobal 
instability  analysis  to  address  the  issue  of  vortex  breakdown  in  this  or  high-angle-of-attack  flow  over 
a  delta- wing  using  this  novel  methodology.  In  the  problem  at  hand,  from  a  physical  point  of  view,  the 
assumptions  of  inviscid  theory  are  increasingly  applicable  as  the  flow  Mach  number  increases.  Although 
subsonic  BiGlobal  calculations  are  the  easiest  to  perform  numerically,  neglecting  all  but  the  streamwise 
velocity  component  in  this  case  is  highly  questionable  and  may  confuse  potential  discrepancies  of 
BiGlobal  instability  analysis  results  with  those  of  experiment /DNS.  Consequently,  further  efforts  on 
BiGlobal  instability  analysis  on  the  elliptic  cone  should  focus  on  the  supersonic/hypersonic  regimes 
alone.  Further  progress  in  this  area  is  expected  to  be  made  by  close  interaction  with  experiment  and/or 
DNS,  so  that  physically  interesting  configurations  and  parameter  ranges  are  isolated  and  studied  in 
detailed  in  the  framework  of  the  present  analysis. 
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Appendix  A.  The  elliptic  confocal  coordinate  system 

Take  the  conformal  transformation 


W=2{Z+z)’ 

where  z  and  w  are  complex  variables.  In  polar  notation, 


which  leads  to 


z  —  relv 


w  —  —  \  re 171  +  -e~ir) 


(Al) 

(A  2) 

(A3) 

(A  4) 
(A  5) 

/ 

r  —  e1^,  (A  6) 

which  results  in  the  relationship  between  the  cartesian  and  the  elliptic  confocal  coordinate  system 


Introducing  w  —  4  (x  +  iy) ,  where  c  is  a  real  parameter,  leads  to 
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sin  T]. 

Key  to  the  transformation  is  the  introduction  of 


x  —  c  cosh  £  cos  rj 
y  —  c  smh  £  sin  r/, 

c  being  identified  as  the  focal  distance  of  the  ellipse. 

Combining  (A3)  and  (A 6)  one  arrives  at 

x  +  iy  —  ccosh(£  +  iry)  (A  9) 

which  may  be  used  to  derive  the  inverse  transformation,  necessary  for  the  numerical  work  in  the 
O^rj  plane  given  the  basic  flow  results  on  the  Oxy  plane, 


(A  7) 
(AS) 


£  +  i  rj  —  —  cosh  1(x  +  iy). 


(A  10) 


The  semi-axes  of  the  ellipse  can  be  obtained  by  introducing  into  (A  7-A  8)  the  limits  y  —  0  and 
V  =  *7% 


a  —  ccosh£,  (A  11) 

&  =  csinh£.  (A  12) 

Rearranging  (A11-A12)  the  well-known  relationship  between  the  focal  distance  and  the  semiaxes 
of  the  ellipse 


2  1 2 
—  a  —  b 


(A  13) 
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may  also  be  obtained.  Also,  combining  (A  11-A  12)  the  equation  of  the  ellipse  fw  —  cnst.  describing 
the  wall  boundary  is  obtained, 


tanh^ 

I  lii'.  may  be  rearranged  to  arrive  at 


sinh  £ u 
cosh  fu 


b  e^w  —  e 


(A  14) 


.  1  .  a  +  b  .  a  +  b 

Zw  =  77  hi - -  =  In - , 

2  a  —  b  c 


(A  15) 


which  is  the  equation  of  the  wall  boundary  in  the  elliptic  confocal  coordinate  system  in  terms  of  the 
semiaxes  and  focal  distance  of  the  ellipse. 
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Appendix  B.  The  temporal  and  spatial  eigenvalue  problems  in  elliptic  confocal 

coordinates 

Given  basic  flow  quantities  and  their  derivatives  on  the  elliptic  confocal  system,  the  eigenvalue 
problems  are  defined  by  the  matrices  T\  —  T12  and  S\  —  S\2-, 


Ti  —  (3  7  ppw 
T2  —  p'yppw 

T3  =  {Ppp^i v  -  P'yp^pw  -  2/3'yppw£  ) 

T4  =  (pppvw  -  P'ypvpw  -  2^7 ppwr^j 

T5  —  ^/33  p2  w3  —  h2/33  r /ppw  ^ 

t6  =  1PP 
t7  -  1PP 

Ts  =  (ppf  -7 P(p) 

T9  =  (ppri  -7  Pnp) 

Tio  =  (3/32p2w  -  h2  /327pp) 

T41  =  —3  [3~p2  w 
Tl2  =  P2 

<5i  =  w  7  PP 

52  —  ioppp 

53  =  w  (pp?  -  7P^p) 

^4  =  w  (ppv  -7 P,,p) 

55  =  w3  p2 

<56  =7  Ppw 
Sr  =7  ppw 

<58  =  (pp£  to  —  7  p^  pw  —  2  7  pp  j 

<59  =  (pPrtW  —  7 p,j pro  -  2-fppw^ 

Sio  —  3  co  2p2  w 

<5n  =  oj  [H2  7pp  —3  p2 

<512  =  (p2  to3  —  ti2  7  pp  w  j 

in  the  temporal  and  spatial  case,  respectively. 
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Appendix  C.  Calculation  of  differential  operators 
in  arbitrary  orthonormal  coordinate  systems 

(a)  Preliminaries 

The  first  step  required  by  differential  form  theory  [11,  3]  is  to  calculate  the  dual  basis  (uis,  u>n,  ujs) 
associated  with  a  vector  space  according  to  the  correspondence 


vector  components 

$ 

(vs,Vn,Vff) 
basis  vectors 

t 

(®si  Gfn 


1- Forms 

t 

(uF,chn,w0) 
basis  1-Forms 

t 

( ds ,  dn ,  d6). 


The  1-Forms  satisfy  the  schematic  relation 


dx% 


r  d 

1 -dxk 

=  £*,  “[eh 

=  4, 


(Cl) 


where  square  brackets  signify  application  of  the  operator  to  the  vector/l-Form  inside  the  bracket 
and  5\  is  the  Kronecker  delta. 

The  second  step,  is  to  use  the  dual  basis  to  express  the  operators  appearing  in  the  equations  of 
motion.  For  example,  the  gradient  operator  is  the  1-Form  df,  defined  as  the  exterior  derivative  of  a 
scalar  field  /, 


df  —  df/dx1  dx\ 


(C  2) 


xl  being  any  chosen  coordinate  system  and  dxl  the  corresponding  basis  1-Forms.  The  definitions 
(C  1)  are  then  used  to  express  the  coefficients  of  the  expansion  in  (C  2)  in  terms  of  the  dual  basis 

jw*  j;  the  coefficients  of  the  result  express  the  gradient  operator  in  an  arbitrary  coordinate  system  x\ 
(i)  Gradient 

For  simplicity  of  exposition  of  the  ideas  of  differential  form  theory  we  calculate  the  gradient  operator 
in  cylindrical  coordinates.  While  this  approach  turns  out  to  be  rather  long-winded  for  this  simple 
problem,  differential  form  theory  can  also  be  applied  in  case  no  analytic  expressions  are  available  to 
link  the  original  and  the  target  coordinate  system;  this  is  precisely  the  case  with  two  of  the  three 
coordinate  systems  discussed  in  §  3(c)  in  the  elliptic  cone. 

In  vector-component  form  the  coordinate  transformation  is  defined  by 

(x  \  /  rcosf  \ 

This  yields 

dx  —  r  cos  f  dr  —  r  sin  f  df. 
dy  —  r  sin  f  dr  +  r  cos  f  df, 
dz  —  dz. 
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The  unknown  1-Forms  (ojr,u)^,uiz)  can  be  expressed  in  terms  of  the  basis  1-Forms  ( dx ,  dy,  dz)  as 

cbr  —  Corx  dx  +  Co'y  dy  +  Corz  dz , 

—  uj%  dx  +  dy  +  uif  dz, 

uz  —  ujz  dx  +  dy  +  uzz  dz. 

where  w*  are  unknown  coefficients.  Application  of  (Cl)  yields  three  3x3  systems  for  the  unknown 

coefficients  w*  of  the  1-Forms. 

For  example, 


( urx  dx  +  ojy  dy  +  u>rz  dz^j  |er  =  1 
[uFx  dx  +  ujy  dy  +  Corz  dz'j  |e<£  =  0, 
[Corx  dx  +  ujy  dy  +  ufz  dz )  |j 


=  0. 


Substituting  the  basis  1-Forms  (dx,dy,dz)  in  terms  of  ( dr,d<f,dz )  and  recalling  (C  1),  i.e.  dr\er 
L,</^[er  =  0, dz^er  —  0,  the  system 


cos  <j)  uorx  +  sin  <f)  ujy  —  1  jr, 

—  sin  (f>  ujx  +  cos  (f>  uiy  —  0, 

«*  =  0 

is  obtained  from  which  ofx  —  1/rcos^,  ujry  —  l/rsmcf,  and  ofz  —  0  follows.  Substitution  of  these 
quantities  into  the  expansion  of  ojr  and  use  of  (C  3-C  5)  completes  the  calculation  of  car.  In  an  analogous 
manner  the  other  components  of  the  dual  basis  are  obtained, 


ur  =  dr 

(C  6) 

lo^  —  r  d(f 

(C  7) 

uz  =  dz. 

(C  8) 

Next,  the  dual  basis  calculated  is  used  to  express  the  gradient  operator.  In  the  chosen  coordinate 
system  the  gradient  operator  takes  the  form 


df  df  df 

df = ir + 4  w + dz 


if  18/ 

dr  r  d< p 


V  +  lfa- 

UZ 


where  the  1-Forms  ( dr,d<f>,dz )  have  been  expressed  in  the  calculated  dual  basis  (C6-C8).  One 
recognizes  that  the  components  of  the  1-Form. 


{dj_  1  d£  d£ 

[  dr1  rd<t>'  dz’ 


express  the  gradient  operator  in  cylindrical  coordinates. 
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(b)  The  gradient  operator  in  the  (eSlenieg)  coordinate  system  of  §  3(c) 

In  the  problem  at  hand  a  vector  can  be  expressed  by  the  unit  vectors  (es,  en .  eg )  and  the  definition 
of  the  coordinate  transformation  (3.23-3.25)  yields 

dx  —  a  cos  <j>  ds  —  (as  +  nbL )  sin  f>  df>  +  bL  cos  <f>  dn , 
dy  —  6  sin  cf  ds  +  ( bs  +  naL )  cos  <f>  d<f>  +  aL  sin  <fi  dn , 
dz  —  L  ds  —  ab  dn , 

and 


^  (a  cos  <fi 

es  =  V£  (  bSL 


—  VS[a  cos  <f>  dx  +  b  sin  <p  dy  +  L  dz] 

—  VS{  | a2  cos2  4>  +  b2  sin2  <f>  +  L2  ds  +  s  \b2  —  a 2  sin  <f>  cos  <f>  d<j> j, 


^  /  Lb  cos  (f  \ 

en  =  -7=  Lasincf)  => 

^  V  ~ab  ) 

uon  —  VN[Lb cos  (f  dx  —  La  sin  cf  dy  —  ab  dz] 

—  y/N  ^nLf2\a2  +  if2  sin  <f>  cos  <p  d<f>  +  L2  |a2  sin2  <j>  +  b2  cos2  <f>  +  1  dnj,  (CIO) 


^  /  a[b2  +  L2]  sin  (f  \ 

e$  —  —j=  —  b[a2  +  L2]  cos  4>  => 

^0  \  L[b2  —  a2]  sin  cfcos  <p  ) 

V  —  \/©|a[62  +  L2]  sin  <f>  dx  —  b[a2  +  L2]  cos  <f>  dy  +  L[b2  —  a2]  sin  <j)  cos  <p  dzj, 

=  —  \/@ja^&2  +  L2^ |as  +  nbL  sin2  </>  + 6  ^a2  +  L2^hs  +  naL  cos2  j  d</>.  (Cll) 

The  3x3  system  for  the  unknown  basis  1-Forms  (ds,dn,  d(/>)  can  be  solved  for  in  terms  of  the  dual 
basis  (ojs,  uin,  Coe).  In  order  to  simplify  the  algebra  we  write 

cbs  —  ujss  ds  +  dcf)  +  0  dn,  (C  12) 

uin  —  0  ds  +  u!1),  d<f  +  c jj™  dn ,  (C  13) 

uis  —  0  ds  +  lo^  d(f>  +  0  dn,  (C  14) 

where 

co')  —  VS  |a2  cos2  (f>  +  b2  sin2  cf  +  L2  , 

=  VSs  ^b2  —  a2  sin </> cos  <j>\ 
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uj, 


uj, 


—  Vn nL2  |a2  +  b 2  sin  (j>  cos  <f>. 

—  Vn L 2 1  a2  sin2  4>  +  b2  cos2  4>  +  1  , 


=  -VQ\a 


(b2  +  L2)  [c 


|as  +  nbL 
The  solution  of  system  (C  12-C  14)  is 


sin 


4>  +  b(a2  +  L2^J  bs 


+  naL 


ds  —  - U!S - , 

wi  “H 


d(f>  =  ~n  UJ  \ 


1  /  )U 
dn  —  — uon - 


Rearranging  the  gradient  operator  (C  2) 

df  —  (df  /ds)  ds  +  (df  /d<f)  d<f  +  (df  /dn)  dn 


(C  15) 
(C  16) 
(C  17) 

(C  18) 


by  substituting  (C  15-C  17)  into  (C  18)  the  gradient  operator  in  the  the  coordinate  system  (s,n,  9) 
can  be  obtained, 


,,  1  df  r  1 

df  —  +  —Q 

Ujf  OS 


UJ 


<t> 


UJ 


<t> 


wfcnj  u 2 


d<f)  uj %  dn 


(C  19) 
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Appendix  D.  Iterative  scheme  for  the  leading  eigenvalue 

At  a  given  streamwise  position  Q  equations  (5. 1-5.2)  are  discretised  at  the  ( N  +  1)  spectral  col¬ 
location  grid  points  yielding  a  linear  inhomogeneous  algebraic  system  which  can  be  expressed  in  the 
form 


AX  =  bo 


A\  A2 

’  <P  ' 

A3  A4 

.  P  . 

=  b, 


(Dl) 


where  <p  —  {<p o,  ip\, ...,  <Pn)T,  P  —  (PthPi-  ■  ■■,Pn)T  ■  The  Nth  and  2JVth  lines  of  matrix  A  are  reserved  for 
the  boundary  conditions.  All  but  two  elements  of  the  vector  b  are  zero.  The  nonzero  elements  appear 
at  the  JVth  and  21Vth  positions  corresponding  to  the  asymptotic  values  of  (p  and  p  at  the  far-held 
position  r]maxi  i.e.,  conditions  (5. 4-5. 5). 

Starting  with  a  suitable  initial  guess  of  the  required  eigenvalue,  the  algebraic  system  (D  1)  is  solved 
iteratively,  by  means  of  a  Newton  iteration  scheme  [17],  until  the  impermeability  condition  on  the 
surface  (i.e.,  <p(0)  =  0)  is  satisfied.  The  above  described  scheme  is  able  to  perform  temporal  as  well 
as  spatial  instability  calculations  giving  spectrally  accurate  results  at  the  same  level  of  computing 
effort.  A  similar  iterative  scheme  was  proposed  by  Malik  [19]  for  the  solution  of  viscous,  compressible 
instability  equations  on  a  real  grid.  Standard  numerical  library  subroutines  were  used  for  the  solution 
of  the  linear  system  (D  1)  and  the  evaluation  of  the  modified  Bessel  functions.  The  non-zero  elements 
of  the  submatrices  A3  and  A4  of  the  system  (D  1)  are 


A\jk  -  Djk  + 


C* 


Wont 


1  +  A  Q  +  Qr)j  \Vq-  —  (3 


-iT0j 


A2jk  7 

_ia2(Wo.  ~P) 

3  jk  j,  1 


2/-2 


1  + 


n2  C 


a2(l  +  A(f  +  CiVj)2 


i (W0j  ~  P) 


7 


A,.,  = 


D 


jk 


jk  (7^)’ 


for  j  —  k.  and 


Atjk  —  Djk, 

A 

= 

for  j  k.  where  j,  k  =  0,  2, ...,  N.  The  elements  of  vector  b  are  given  by 

bj  =  0,  if  j  7^  N  and  j  7^  2 IV, 

bN  =  ^ \Kn+i{f)i )  + 
h  =  M^iaryjl  -  p)Kn{r)i) 

2N  *  [l-M2o(l-/32)]l/2  ’ 


where 


and  j,  k  —  0, 2,  ...,21V. 


Vi  =  ±a[l  -  M£(  1  -  /32)]1/2(1/C *  +  +  Vmax), 


(D2) 

(D3) 

(D4) 

(D5) 

(D6) 

(D7) 

(D8) 

(D9) 

(D  10) 
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N 

O 

II 

iM 

'O 

52  =  0.05 

T— 1 

O 

II 

'O 

52  =  0.15 

C 

c 

C 

c 

48 

0.856328  +i  0.000000 

0.856333  +i  0.000418 

0.856129  +i  0.000355 

0.855075  -i  0.000462 

64 

0.856543  +i  0.000000 

0.856357  +i  0.000407 

0.856382  +i  0.000395 

0.856723  -i  0.000134 

80 

0.856937  +i  0.000000 

0.856356  +i  0.000408 

0.856354  +i  0.000412 

0.856228  +i  0.000550 

96 

0.857174  +i  0.001755 

0.856356  +i  0.000408 

0.856356  +i  0.000408 

0.856324  +i  0.000356 

Table  1.  Dependence  of  the  eigenvalue  c  on  the  number  of  collocation  nodes  and  the  complex  mapping  parameter  62  in  the  case 
of  planar,  compressible,  adiabatic  flow  for  a  =  0.365  (mode  II)  and  Moo  =  3.8. 
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N 

f]max 

l 

Si 

Mode  I 

lo  —  0.0386,  n  —  3 

Mode  II 

cj  =  0.334,  n  —  0 

32 

25 

25 

-0.05 

a 

0.050338  -i  0.003137 

a 

0.367236  -i  0.004590 

64 

25 

25 

-0.05 

0.050338  -i  0.003137 

0.366888  -i  0.004487 

80 

25 

25 

-0.05 

0.050339  -i  0.003136 

0.366903  -i  0.004485 

96 

25 

25 

-0.05 

0.050339  -i  0.003136 

0.366903  -i  0.004485 

80 

40 

40 

-0.04 

0.050339  -i  0.003136 

0.366905  -i  0.004483 

90 

60 

60 

-0.04 

0.050339  -i  0.003136 

0.366902  -i  0.004483 

Table  2.  Dependence  of  the  eigenvalues  on  the  number  of  collocation  nodes  and  the  mappings  parameters,  for  the  case  of 

adiabatic  cone  flow  with  Moo  =  3.8  and  A  =  1,  at  C  =  0.05. 
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Nr] 

fn(x  -  a,y  —  0) 

fs(x  =  a,y  —  0) 

fA(x  -  a,y  -  0) 

10 

10 

-0.901970 

-1.360926 

-0.863940 

20 

20 

-0.901970 

-1.367825 

-0.867146 

40 

40 

-0.901970 

-1.375666 

-0.868566 

80 

80 

-0.901970 

-1.375823 

-0.868586 

Table  3.  Convergence  history  of  numerical  solution  of  the  Poisson  problems  H,  S  and  A,  subject  to  homogeneous  Neumann 
boundary  conditions.  N£  and  Nr]  respectively  denote  the  number  of  collocation  points  along  the  £  and  r]  coordinate  directions. 
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Figure  1.  The  confocal  elliptical  coordinate  system  O^T]. 
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Figure  4.  An  arbitrary  location  #(</>). 
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n 


Figure  5.  One  orthogonal  body-fitted  coordinate  system  which  defines  the  On8  plane  on  which  a  global  instability  analysis  can 

be  performed. 
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Frequencies  cr  and  growth  rates  c\  of  mode  I  (dashed)  and  II  (dash-dotted)  instabilities  as  function  of  the  wavenumber 
=  3.8  boundary-layer  flow  over  an  insulated  flat  plate.  Superimposed,  denoted  by  open  symbols,  are  the  results  of  [8]. 
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Figure  7.  Temporal  instability  results.  Distribution  of  the  temporal  growth  rate  ac»  with  a  for  modes  I  and  II,  C  =  0.05  and 
various  azimuthal  wavenumbers.  For  comparison,  the  results  of  [9]  are  also  shown  indicated  by  symbols. 
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Figure  8.  Calculation  of  geometric  parameters  on  the  instability  analysis  plane 
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Figure  9.  The  structure  of  the  matrix  discretizing  the  spatial  eigenvalue  problem 
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Figure  11.  Numerical  solution  of  the  Poisson  problem  (6.1)  using  (6.3).  Upper:  /(£  =  .  rj)  =  0;  lower:  /^(£  =  £w,ri)  =  0. 
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Figure  12.  Numerical  solution  of  the  Poisson  problem  (6.1)  using  (6.4).  Upper:  /(£  =  .  rj)  =  0;  lower:  /^(£  =  £w,v)  =  0. 
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Figure  13.  Numerical  solution  of  the  Poisson  problem  (6.1)  using  (6.5).  Upper:  /(£  =  .  rj)  =  0;  lower:  /^(£  =  £w,v)  =  0. 
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Figure  17.  A  stationary  mode  at  M  =  0.5.  Upper  p ,  lower  left  p lower  right  pv. 
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Figure  18.  A  stationary  mode  at  M  =  0.5.  Upper  p ,  lower  left  p lower  right  pv. 
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Figure  19.  A  stationary  mode  at  M  =  0.5.  Upper  p ,  lower  left  p lower  right  pv. 
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Figure  21.  A  travelling  mode  -  member  of  family  A  -  at  M  =  4.0.  Upper  p ,  lower  left  p lower  right  pv- 
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Figure  22.  A  travelling  mode  -  member  of  family  B  -  at  M  =  4.0.  Upper  p,  lower  left  p^,  lower  right  pv. 
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Figure  23.  A  travelling  neutral  mode  -  member  of  family  C  -  at  M  =  4.0.  Upper  p,  lower  left  p lower  right  pv. 
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Figure  24.  A  stationary  mode  -  member  of  family  D  -  at  M  =  4.0.  Upper  p ,  lower  left  p lower  right  pv. 
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Figure  25.  A  stationary  neutral  pressure  eigenmode  at  M  =  8.0. 
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Spatial  BiGlobal  Eigenspectrum  at  M  =  8 


Figure  26.  Eigenvalue  spectra  at  M  =  8  and  resolutions  102,202  and  25; 
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Figure  27.  A  stationary  neutral  mode  at  M  =  8.0.  Upper  p,  lower  left  p^,  lower  right  pv. 
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Figure  28.  Disturbance  pressure  travelling  eigenmode  at  M  =  8.0.  Upper,  full  calculation  domain;  lower,  detail  in  the 

neighbourhood  of  the  elliptic  cone. 
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Figure  29.  Disturbance  pressure  travelling  eigenmode  at  M  =  8.0.  Upper,  full  calculation  domain;  lower,  detail  in  the 

neighbourhood  of  the  elliptic  cone. 
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